!@sum DIAG ModelE diagnostic calculations
!@auth G. Schmidt/J. Lerner/R. Ruedy/M. Kelley
!@ver 1.0
C**** AJ(J,N) (ZONAL SUM OVER LONGITUDE AND TIME)
C**** See j_defs for contents
C**** IDACC
C**** CONTENTS OF APJ(J,N) (SUM OVER LONGITUDE AND TIME OF)
C**** 1 P (100 PA) 4 DA
C**** 2 4*P4I (100 PA) (UV GRID) 4 DA
C****
C**** CONTENTS OF AJL(J,L,N) (SUM OVER LONGITUDE AND TIME OF)
C**** See jl_defs for contents
C****
C**** CONTENTS OF ASJL(J,L,N) (SUM OVER LONGITUDE AND TIME OF)
C**** See jls_defs for contents
C****
C**** CONTENTS OF AIJ(I,J,N) (SUM OVER TIME OF)
C**** See ij_defs for contents
C****
C**** CONTENTS OF AIL(I,L,N) (SUM OVER TIME OF)
C**** See il_defs for contents
C****
C**** CONTENTS OF IDACC(N), NUMBER OF ACCUMULATION TIMES OF
C**** 1 SOURCE TERMS (dt: DTSRC)
C**** 2 RADIATION SOURCE TERMS (dt: NRAD*DTsrc)
C**** 3 SURFACE INTERACTION SOURCE TERMS (dt: NDASF*DTsrc+DTsurf)
C**** 4 QUANTITIES IN DIAGA (dt: NDAA*DTsrc+2*DTdyn)
C**** 5 ENERGY NUMBERS IN DIAG4 (DEYERMINED BY NDA4)
C**** 6 KINETIC ENERGY IN DIAG5 FROM DYN'CS (dt: NDA5K*DTsrc+2*DTdyn)
C**** 7 ENERGY IN DIAG5 FROM DYNAMICS (dt: NDA5D*DTsrc)
C**** 8 ENERGY IN DIAG5 FROM SOURCES (DETERMINED BY NDA5S)
C**** 9 WAVE ENERGY IN DIAG7 (dt: 12 HOURS)
C**** 10 ENERGY IN DIAG5 FROM FILTER (DT: NFILTR*DTsrc)
C**** 11 NOT USED
C**** 12 ALWAYS =1 (UNLESS SEVERAL RESTART FILES WERE ACCUMULATED)
C****
MODULE DIAG_LOC 7,1
!@sum DIAG_LOC is a local module for some saved diagnostic calculations
!@auth Gavin Schmidt
USE MODEL_COM
, only : im,imh,jm,lm
IMPLICIT NONE
SAVE
C**** Variables passed from DIAGA to DIAGB
!@var W,TX vertical velocity and in-situ temperature calculations
REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: W
REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: TX
!@var TJL0 zonal mean temperatures prior to advection
REAL*8, ALLOCATABLE, DIMENSION(:,:) :: TJL0
C**** Variables used in DIAG5 calculations
!@var FCUVA,FCUVB fourier coefficients for velocities
REAL*8, ALLOCATABLE, DIMENSION(:,:,:,:) :: FCUVA,FCUVB
C**** Some local constants
!@var JET, LDEX model levels for various pressures
!@var LUPA,LDNA shorthand for above/below levels
!@var PMO,PLO,PM,PL some shorthand pressure level
INTEGER :: JET
INTEGER, DIMENSION(3) :: LDEX
REAL*8, DIMENSION(LM) :: LUPA,LDNA
REAL*8, DIMENSION(LM) :: PMO,PLO
REAL*8, DIMENSION(LM+1) :: PM,PL
END MODULE DIAG_LOC
SUBROUTINE ALLOC_DIAG_LOC(grid) 1,5
USE DOMAIN_DECOMP, only : GET
USE DOMAIN_DECOMP, only : DYN_GRID
USE MODEL_COM, only : im,imh,lm
USE DIAG_LOC, only : W,TX,TJL0,FCUVA,FCUVB
IMPLICIT NONE
LOGICAL, SAVE :: init=.false.
INTEGER :: J_1H , J_0H
INTEGER :: IER
TYPE(DYN_GRID) :: grid
If (init) Then
Return ! Only invoke once
End If
init = .true.
CALL GET(grid, J_STRT_HALO=J_0H, J_STOP_HALO=J_1H)
ALLOCATE( W(IM, J_0H:J_1H, LM), TX(IM, J_0H:J_1H, LM),
& STAT = IER)
ALLOCATE( TJL0(J_0H:J_1H, LM),
& STAT = IER)
ALLOCATE( FCUVA(0:IMH, J_0H:J_1H, LM, 2),
& FCUVB(0:IMH, J_0H:J_1H, LM, 2),
& STAT = IER)
RETURN
END SUBROUTINE ALLOC_DIAG_LOC
SUBROUTINE DIAGA 1,35
!@sum DIAGA accumulate various diagnostics during dynamics
!@auth Original Development Team
!@ver 1.0
USE CONSTANT, only : grav,rgas,kapa,lhe,sha,bygrav,bbyg,gbyrb,tf
* ,rvap,gamd,teeny,undef
USE MODEL_COM, only : im,imh,fim,byim,jm,jeq,lm,ls1,idacc,ptop
* ,pmtop,psfmpt,mdyn,mdiag,sig,sige,dsig,zatmo,WM,ntype,ftype
* ,u,v,t,p,q
USE GEOM, only : areag,cosp,dlat,dxv,dxyn,dxyp,dxys,dxyv,dyp,fcor
* ,imaxj,ravpn,ravps,sinp,bydxyv
USE RADNCB, only : rqt,lm_req
USE DAGCOM, only : aj,areg,jreg,apj,ajl,asjl,ail,j50n,j70n,j5nuv
* ,j5suv,j5s,j5n,aij,ij_dtdp,ij_dsev,ij_phi1k,ij_pres,ij_puq
* ,ij_pvq,ij_slp,ij_t850,ij_t500,ij_t300,ij_q850,ij_q500
* ,ij_RH1,ij_RH850,ij_RH500,ij_RH300
* ,ij_q300,ij_ujet,ij_vjet,j_tx1,j_tx,j_qp,j_dtdjt,j_dtdjs
* ,j_dtdgtr,j_dtsgst,j_rictr,j_rostr,j_ltro,j_ricst,j_rosst
* ,j_lstr,j_gamm,j_gam,j_gamc,lstr,il_ueq,il_veq,il_weq,il_teq
* ,il_qeq,il_w50n,il_t50n,il_u50n,il_w70n,il_t70n,il_u70n
* ,kgz_max,pmb,ght,jl_dtdyn,jl_zmfntmom,jl_totntmom,jl_ape
* ,jl_uepac,jl_vepac,jl_uwpac,jl_vwpac,jl_wepac,jl_wwpac
* ,jl_epflxn,jl_epflxv,ij_p850,z_inst,rh_inst,t_inst,plm
USE DYNAMICS, only : pk,phi,pmid,plij, pit,sd,pedn
USE PBLCOM, only : tsavg
USE DIAG_LOC, only : w,tx,lupa,ldna,jet,tjl0
USE DOMAIN_DECOMP, only : GET, CHECKSUM, HALO_UPDATE, GRID
USE DOMAIN_DECOMP, only : SOUTH, NORTH
IMPLICIT NONE
REAL*8, DIMENSION(LM) :: GMEAN
REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO) ::
& TIL,UI,UMAX,PI,EL,RI,DUDVSQ
REAL*8, DIMENSION(NTYPE,GRID%J_STRT_HALO:GRID%J_STOP_HALO) ::
& SPTYPE
REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM) ::
& THJL,THSQJL,SPI,PHIPI,TPI
REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM-1) ::
& SDMEAN
REAL*8, DIMENSION(IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO) ::
& PUV
REAL*8, DIMENSION(LM_REQ) :: TRI
REAL*8, DIMENSION(IM) :: THSEC,PSEC,SQRTP,PDA
CHARACTER*16 TITLE
REAL*8, PARAMETER :: ONE=1.,P1000=1000.
INTEGER :: I,IM1,J,K,L,JR,LDN,LUP,
& IP1,LM1,LP1,LR,MBEGIN,IT
REAL*8 THBAR ! external
REAL*8 ::
& BBYGV,BDTDL,BYSDSG,CDTDL,DLNP,DLNP12,DLNP23,DBYSD,
& DLNS,DP,DS,DT2,DTHDP,DU,DUDP,DUDX,DV,DXYPJ,ELX,
* ESEPS,FPHI,GAMC,GAMM,GAMX,GMEANL,P1,P4,P4I,
& PDN,PE,PEQ,PEQM1,PEQM2,PHIRI,PIBYIM,PIJ,PITIJ,PITMN,
* PKE,PL,PRT,PU4I,PUP,PUV4I,PV4I,PVTHP,
* QLH,ROSSX,SDMN,SDPU,SMALL,SP,SP2,SS,T4,THETA,THGM,THMN,TPIL,
* TZL,UAMAX,UMN,UPE,VPE,X,Z4,THI,TIJK,QIJK
LOGICAL qpress,qabove
INTEGER nT,nQ,nRH
REAL*8, PARAMETER :: EPSLON=1.
INTEGER, PARAMETER ::
* I150E = IM*(180+150)/360+1, ! WEST EDGE OF 150 EAST
* I110W = IM*(180-110)/360+1, ! WEST EDGE OF 110 WEST
* I135W = IM*(180-135)/360+1 ! WEST EDGE OF 135 WEST
REAL*8 QSAT
INTEGER :: J_0, J_1, J_0S, J_1S, J_0STG, J_1STG
CALL GETTIME(MBEGIN)
CALL GET(grid, J_STRT=J_0, J_STOP=J_1,
& J_STRT_SKP=J_0S, J_STOP_SKP=J_1S,
& J_STRT_STGR=J_0STG, J_STOP_STGR=J_1STG)
IDACC(4)=IDACC(4)+1
BYSDSG=1./(1.-SIGE(LM+1))
DLNP12=LOG(.75/.35)
DLNP23=LOG(.35/.1)
C****
C**** FILL IN HUMIDITY AND SIGMA DOT ARRAYS AT THE POLES
C****
IF(GRID%HAVE_SOUTH_POLE) THEN
DO L=1,LM
DO I=2,IM
Q(I,1,L)=Q(1,1,L)
END DO
END DO
ENDIF ! GRID%HAVE_SOUTH_POLE
IF(GRID%HAVE_NORTH_POLE) THEN
DO L=1,LM
DO I=2,IM
Q(I,JM,L)=Q(1,JM,L)
END DO
END DO
ENDIF ! GRID%HAVE_NORTH_POLE
C****
C**** CALCULATE PK AND TX, THE REAL TEMPERATURE
C****
IF(GRID%HAVE_SOUTH_POLE) THEN
DO L=1,LM
TX(1,1,L)=T(1,1,L)*PK(L,1,1)
DO I=2,IM
T(I,1,L)=T(1,1,L)
TX(I,1,L)=TX(1,1,L)
END DO
END DO
ENDIF ! GRID%HAVE_SOUTH_POLE
IF(GRID%HAVE_NORTH_POLE) THEN
DO L=1,LM
TX(1,JM,L)=T(1,JM,L)*PK(L,1,JM)
DO I=2,IM
T(I,JM,L)=T(1,JM,L)
TX(I,JM,L)=TX(1,JM,L)
END DO
END DO
ENDIF ! GRID%HAVE_NORTH_POLE
DO L=1,LM
DO J=J_0S,J_1S
DO I=1,IM
TX(I,J,L)=T(I,J,L)*PK(L,I,J)
END DO
END DO
END DO
C****
C**** CALCULATE PUV, THE MASS WEIGHTED PRESSURE
C****
CALL CHECKSUM(grid, P, 224, "DIAG.f")
CALL HALO_UPDATE(grid, P, FROM=SOUTH)
DO J=J_0STG,J_1STG
I=IM
DO IP1=1,IM
PUV(I,J)=RAVPN(J-1)*( * RAVPS( J)*( I=IP1
END DO
END DO
C****
C**** J LOOPS FOR ALL PRIMARY GRID ROWS
C****
DO J=J_0,J_1
DXYPJ=DXYP(J)
C**** NUMBERS ACCUMULATED FOR A SINGLE LEVEL
PI(J)=0.
SPTYPE(:,J)=0
DO I=1,IMAXJ(J)
JR=JREG(I,J)
DO IT=1,NTYPE
SPTYPE(IT,J)=SPTYPE(IT,J)+FTYPE(IT,I,J)
AJ(J,J_TX1,IT)=AJ(J,J_TX1,IT)+(TX(I,J,1)-TF)*FTYPE(IT,I,J)
END DO
AREG(JR,J_TX1)=AREG(JR,J_TX1)+(TX(I,J,1)-TF)*DXYPJ
PI(J)=PI(J)+P(I,J)
AIJ(I,J,IJ_PRES)=AIJ(I,J,IJ_PRES)+ P(I,J)
AIJ(I,J,IJ_SLP)=AIJ(I,J,IJ_SLP)+((P(I,J)+PTOP)*(1.+BBYG
* *ZATMO(I,J)/TSAVG(I,J))**GBYRB-P1000)
AIJ(I,J,IJ_RH1)=AIJ(I,J,IJ_RH1)+Q(I,J,1)/QSAT(TX(I,J,1),LHE,
* PMID(1,I,J))
END DO
APJ(J,1)=APJ(J,1)+PI(J)
C**** CALCULATE GEOPOTENTIAL HEIGHTS AT SPECIFIC MILLIBAR LEVELS
DO I=1,IMAXJ(J)
K=1
L=1
rh_inst(:,i,j) = undef ; t_inst(:,i,j) = undef
z_inst(:,i,j) = undef
172 L=L+1
PDN=PMID(L-1,I,J)
PL=PMID(L,I,J)
IF (PMB(K).LT.PL.AND.L.LT.LM) GO TO 172
C**** Select pressure levels on which to save temperature and humidity
C**** Use masking for 850 mb temp/humidity
174 qpress = .false.
qabove = pmb(k).le.pedn(l-1,i,j)
SELECT CASE (NINT(PMB(K)))
CASE (850) ! 850 mb
nT = IJ_T850 ; nQ = IJ_Q850 ; nRH = IJ_RH850 ; qpress=.true.
if (.not. qabove) qpress = .false.
if (qpress) aij(i,j,ij_p850) = aij(i,j,ij_p850) + 1.
CASE (500) ! 500 mb
nT = IJ_T500 ; nQ = IJ_Q500 ; nRH = IJ_RH500 ; qpress=.true.
CASE (300) ! 300 mb
nT = IJ_T300 ; nQ = IJ_Q300 ; nRH = IJ_RH300 ; qpress=.true.
END SELECT
C**** calculate geopotential heights + temperatures
IF (ABS(TX(I,J,L)-TX(I,J,L-1)).GE.EPSLON) THEN
BBYGV=(TX(I,J,L-1)-TX(I,J,L))/(PHI(I,J,L)-PHI(I,J,L-1))
AIJ(I,J,IJ_PHI1K-1+K)=AIJ(I,J,IJ_PHI1K-1+K)+(PHI(I,J,L)
* -TX(I,J,L)*((PMB(K)/PL)**(RGAS*BBYGV)-1.)/BBYGV-GHT(K)
* *GRAV)
IF (qabove) then
TIJK=(TX(I,J,L)-TF
* +(TX(I,J,L-1)-TX(I,J,L))*LOG(PMB(K)/PL)/LOG(PDN/PL))
Z_inst(K,I,J)=(PHI(I,J,L)
* -TX(I,J,L)*((PMB(K)/PL)**(RGAS*BBYGV)-1.)/BBYGV-GHT(K)
* *GRAV)
END IF
ELSE
AIJ(I,J,IJ_PHI1K-1+K)=AIJ(I,J,IJ_PHI1K-1+K)+(PHI(I,J,L)
* -RGAS*TX(I,J,L)*LOG(PMB(K)/PL)-GHT(K)*GRAV)
IF (qabove) then
TIJK=TX(I,J,L)-TF
Z_inst(K,I,J)=(PHI(I,J,L)
* -RGAS*TX(I,J,L)*LOG(PMB(K)/PL)-GHT(K)*GRAV)
END IF
END IF
if (qabove) then
QIJK=Q(I,J,L)+( RH_inst(K,I,J)=QIJK/qsat(TIJK+TF,LHE,PMB(K))
T_inst(K,I,J) =TIJK
if (qpress) then
AIJ(I,J,nT)=AIJ(I,J,nT)+TIJK
AIJ(I,J,nQ)=AIJ(I,J,nQ)+QIJK
AIJ(I,J,nRH)=AIJ(I,J,nRH)+QIJK/qsat(TIJK+TF,LHE,PMB(K))
end if
end if
C****
IF (K.LT.KGZ_max) THEN
K=K+1
IF (PMB(K).LT.PL.AND.L.LT.LM) GO TO 172
GO TO 174
END IF
END DO
END DO
C**** ACCUMULATION OF TEMP., POTENTIAL TEMP., Q, AND RH
DO J=J_0,J_1
DXYPJ=DXYP(J)
DO L=1,LM
TPI(J,L)=0.
PHIPI(J,L)=0.
SPI(J,L)=0.
THI=0.
DBYSD=DSIG(L)*BYSDSG
DO I=1,IMAXJ(J)
JR=JREG(I,J)
PIJ=PLIJ(L,I,J)
DO IT=1,NTYPE
AJ(J,J_TX,IT)=AJ(J,J_TX,IT)+(TX(I,J,L)-TF)*FTYPE(IT,I,J)
* *DBYSD
AJ(J,J_QP,IT)=AJ(J,J_QP,IT)+( * *DSIG(L)*FTYPE(IT,I,J)
END DO
AREG(JR,J_QP)=AREG(JR,J_QP)+( * *DXYPJ
AREG(JR,J_TX)=AREG(JR,J_TX)+(TX(I,J,L)-TF)*DBYSD*DXYPJ
TPI(J,L)=TPI(J,L)+(TX(I,J,L)-TF)*PIJ
PHIPI(J,L)=PHIPI(J,L)+PHI(I,J,L)*PIJ
SPI(J,L)=SPI(J,L)+T(I,J,L)*PIJ
THI=THI+T(I,J,L)
END DO
AJL(J,L,JL_DTDYN)=AJL(J,L,JL_DTDYN)+THI-TJL0(J,L)
END DO
END DO
C****
C**** NORTHWARD GRADIENT OF TEMPERATURE: TROPOSPHERIC AND STRATOSPHERIC
C****
CALL CHECKSUM(grid, TX, 354, "DIAG.f")
CALL HALO_UPDATE(grid, TX, FROM=SOUTH)
CALL HALO_UPDATE(grid, TX, FROM=NORTH)
DO J=J_0S,J_1S
C**** MEAN TROPOSPHERIC NORTHWARD TEMPERATURE GRADIENT
DO L=1,LS1-1
DO I=1,IM
DO IT=1,NTYPE
AJ(J,J_DTDJT,IT)=AJ(J,J_DTDJT,IT)+(TX(I,J+1,L)-TX(I,J-1,L))
* *FTYPE(IT,I,J)*DSIG(L)
END DO
END DO
END DO
C**** MEAN STRATOSPHERIC NORTHWARD TEMPERATURE GRADIENT
DO L=LS1,LSTR
DO I=1,IM
DO IT=1,NTYPE
AJ(J,J_DTDJS,IT)=AJ(J,J_DTDJS,IT)+(TX(I,J+1,L)-TX(I,J-1,L))
* *FTYPE(IT,I,J)*DSIG(L)
END DO
END DO
END DO
END DO
C****
C**** STATIC STABILITIES: TROPOSPHERIC AND STRATOSPHERIC
C****
DO J=J_0,J_1
DXYPJ=DXYP(J)
C**** OLD TROPOSPHERIC STATIC STABILITY
DO I=1,IMAXJ(J)
JR=JREG(I,J)
SS=( DO IT=1,NTYPE
AJ(J,J_DTDGTR,IT)=AJ(J,J_DTDGTR,IT)+SS*FTYPE(IT,I,J)
END DO
AREG(JR,J_DTDGTR)=AREG(JR,J_DTDGTR)+SS*DXYPJ
AIJ(I,J,IJ_DTDP)=AIJ(I,J,IJ_DTDP)+SS
END DO
C**** OLD STRATOSPHERIC STATIC STABILITY (USE LSTR as approx 10mb)
DO I=1,IMAXJ(J)
JR=JREG(I,J)
SS=( * +teeny)
DO IT=1,NTYPE
AJ(J,J_DTSGST,IT)=AJ(J,J_DTSGST,IT)+SS*FTYPE(IT,I,J)
END DO
AREG(JR,J_DTSGST)=AREG(JR,J_DTSGST)+SS*DXYPJ
END DO
C****
C**** NUMBERS ACCUMULATED FOR THE RADIATION EQUILIBRIUM LAYERS
C****
DO LR=1,LM_REQ
TRI(LR)=0.
DO I=1,IMAXJ(J)
TRI(LR)=TRI(LR)+RQT(LR,I,J)
END DO
ASJL(J,LR,1)=ASJL(J,LR,1)+(TRI(LR)-TF*IMAXJ(J))
END DO
PHIRI=0.
DO I=1,IMAXJ(J)
PHIRI=PHIRI+(PHI(I,J,LM)+RGAS*.5*(TX(I,J,LM)+RQT(1,I,J))
* *LOG((SIG(LM)*PSFMPT+PTOP)/PLM(LM+1)))
END DO
ASJL(J,1,2)=ASJL(J,1,2)+PHIRI
PHIRI=PHIRI+RGAS*.5*(TRI(1)+TRI(2))*DLNP12
ASJL(J,2,2)=ASJL(J,2,2)+PHIRI
PHIRI=PHIRI+RGAS*.5*(TRI(2)+TRI(3))*DLNP23
ASJL(J,3,2)=ASJL(J,3,2)+PHIRI
END DO
C****
C**** RICHARDSON NUMBER , ROSSBY NUMBER , RADIUS OF DEFORMATION
C****
C**** NUMBERS ACCUMULATED OVER THE TROPOSPHERE
DO J=J_0STG,J_1STG
DUDVSQ(J)=0.
UMAX(J)=0.
DO I=1,IM
DU=U(I,J,LS1-1)-U(I,J,1)
DV=V(I,J,LS1-1)-V(I,J,1)
DUDVSQ(J)=DUDVSQ(J)+(DU*DU+DV*DV)*PUV(I,J)
END DO
END DO
CALL CHECKSUM(grid, DUDVSQ, 438, "DIAG.f")
CALL HALO_UPDATE(grid, DUDVSQ, FROM=NORTH)
DO J=J_0S,J_1S
PIBYIM=PI(J)*BYIM
DLNP=LOG((SIG(1)*PIBYIM+PTOP)/(SIG(LS1-1)*PIBYIM+PTOP))
DLNS=LOG(SPI(J,LS1-1)/SPI(J,1))
DS=SPI(J,LS1-1)-SPI(J,1)
EL(J)=SQRT(DLNS/DLNP)
RI(J)=DS*DLNP/(.5*(DUDVSQ(J)+DUDVSQ(J+1)))
END DO
DO L=1,LS1-1
DO J=J_0STG,J_1STG
UI(J)=0.
DO I=1,IM
UI(J)=UI(J)+U(I,J,L)
END DO
END DO
CALL CHECKSUM(grid, UI, 457, "DIAG.f")
CALL HALO_UPDATE(grid, UI, FROM=NORTH)
DO J=J_0S,J_1S
UAMAX=ABS(UI(J)+UI(J+1))
IF (UAMAX.GT.UMAX(J)) UMAX(J)=UAMAX
END DO
END DO
DO J=J_0S,J_1S
ROSSX=DYP(J)/(DXYP(J)*SINP(J))
ELX=1./SINP(J)
DO IT=1,NTYPE
AJ(J,J_RICTR,IT)=AJ(J,J_RICTR,IT)+RI(J) *SPTYPE(IT,J)
AJ(J,J_ROSTR,IT)=AJ(J,J_ROSTR,IT)+UMAX(J)*SPTYPE(IT,J)*ROSSX
AJ(J,J_LTRO ,IT)=AJ(J,J_LTRO ,IT)+EL(J) *SPTYPE(IT,J)*ELX
END DO
END DO
C**** NUMBERS ACCUMULATED OVER THE LOWER STRATOSPHERE
C**** LSTR is approx. 10mb level. This maintains consistency over
C**** the different model tops
DO J=J_0STG,J_1STG
DUDVSQ(J)=0.
UMAX(J)=0.
DO I=1,IM
DU=U(I,J,LSTR)-U(I,J,LS1-1)
DV=V(I,J,LSTR)-V(I,J,LS1-1)
DUDVSQ(J)=DUDVSQ(J)+(DU*DU+DV*DV)*PUV(I,J)
END DO
END DO
CALL CHECKSUM(grid, DUDVSQ, 487, "DIAG.f")
CALL HALO_UPDATE(grid, DUDVSQ, FROM=NORTH)
DO J=J_0S,J_1S
PIBYIM=PI(J)*BYIM
DLNP=LOG((SIG(LS1-1)*PIBYIM+PTOP)/(SIG(LSTR)*PSFMPT+PTOP))
DLNS=LOG(SPI(J,LSTR)/SPI(J,LS1-1))
DS=SPI(J,LSTR)-SPI(J,LS1-1)
EL(J)=SQRT(DLNS/DLNP)
RI(J)=DS*DLNP/(.5*(DUDVSQ(J)+DUDVSQ(J+1)))
END DO
DO L=LS1,LSTR
DO J=J_0STG,J_1STG
UI(J)=0.
DO I=1,IM
UI(J)=UI(J)+U(I,J,L)
END DO
END DO
CALL CHECKSUM(grid, UI, 506, "DIAG.f")
CALL HALO_UPDATE(grid, UI, FROM=NORTH)
DO J=J_0S,J_1S
UAMAX=ABS(UI(J)+UI(J+1))
IF (UAMAX.GT.UMAX(J)) UMAX(J)=UAMAX
END DO
END DO
DO J=J_0S,J_1S
ROSSX=DYP(J)/(DXYP(J)*SINP(J))
ELX=1./SINP(J)
DO IT=1,NTYPE
AJ(J,J_RICST,IT)=AJ(J,J_RICST,IT)+RI(J) *SPTYPE(IT,J)
AJ(J,J_ROSST,IT)=AJ(J,J_ROSST,IT)+UMAX(J)*SPTYPE(IT,J)*ROSSX
AJ(J,J_LSTR ,IT)=AJ(J,J_LSTR ,IT)+EL(J) *SPTYPE(IT,J)*ELX
END DO
END DO
C****
C**** MEAN TROPOSPHERIC LAPSE RATES: MOIST CONVECTIVE, ACTUAL,
C**** DRY ADIABATIC
C****
X=RGAS*LHE*LHE/(SHA*RVAP)
DO J=J_0,J_1
GAMM=0.
DO L=1,LS1-1
TZL=TPI(J,L)/PI(J)+TF
PRT=(SIG(L)*PI(J)*BYIM+PTOP)*RGAS*TZL
ESEPS=QSAT(TZL,LHE,ONE)
GAMM=GAMM+(PRT+LHE*ESEPS)/(PRT+X*ESEPS/TZL)*DSIG(L)
END DO
GAMX=(TPI(J,1)-TPI(J,LS1-1))/(PHIPI(J,LS1-1)-PHIPI(J,1))
DO IT=1,NTYPE
AJ(J,J_GAMM,IT)=AJ(J,J_GAMM,IT)+GAMM*SPTYPE(IT,J)
AJ(J,J_GAM ,IT)=AJ(J,J_GAM ,IT)+GAMX*SPTYPE(IT,J)
END DO
END DO
C**** DRY ADIABATIC LAPSE RATE
DO J=J_0,J_1
TPIL=0.
DO L=1,LS1-1
TPIL=TPIL+TPI(J,L)*DSIG(L)
END DO
TIL(J)=TPIL/(PI(J)*(SIGE(1)-SIGE(LS1)))
END DO
CALL CHECKSUM(grid, TIL, 551, "DIAG.f")
CALL HALO_UPDATE(grid, TIL, FROM=NORTH)
CALL HALO_UPDATE(grid, TIL, FROM=SOUTH)
DO J=J_0S,J_1S
X=SINP(J)*GRAV/(COSP(J)*RGAS*2.*DLAT)
DT2=TIL(J+1)-TIL(J-1)
GAMC=GAMD+X*DT2/(TIL(J)+TF)
DO IT=1,NTYPE
AJ(J,J_GAMC,IT)=AJ(J,J_GAMC,IT)+GAMC*SPTYPE(IT,J)
END DO
END DO
C****
C**** EASTWARD TRANSPORTS
C****
CALL CHECKSUM(grid, U, 567, "DIAG.f")
CALL HALO_UPDATE(grid, U, FROM=NORTH)
DO L=1,LM
DO J=J_0S,J_1S
I=IM
DO IP1=1,IM
AIJ(I,J,IJ_PUQ)=AIJ(I,J,IJ_PUQ)+(PLIJ(L,I,J)+PLIJ(L,IP1,J))*
* ( I=IP1
END DO
END DO
END DO
C****
C**** MOMENTUM, KINETIC ENERGY, NORTHWARD TRANSPORTS, ANGULAR MOMENTUM
C****
!Not necessary here, done above CALL CHECKSUM(grid, P, 584, "DIAG.f")
!Not necessary here, done above CALL HALO_UPDATE(grid, P, FROM=SOUTH)
CALL CHECKSUM(grid, PLIJ, 586, "DIAG.f")
CALL HALO_UPDATE(grid, PLIJ, FROM=SOUTH)
!Not necessary here, done above CALL CHECKSUM(grid, TX, 588, "DIAG.f")
!Not necessary here, done above CALL HALO_UPDATE(grid, TX, FROM=SOUTH)
CALL CHECKSUM(grid, PHI, 590, "DIAG.f")
CALL HALO_UPDATE(grid, PHI, FROM=SOUTH)
DO J=J_0STG,J_1STG
P4I=0.
I=IM
DO IP1=1,IM
P4=P(I,J-1)+P(IP1,J-1)+P(I,J)+P(IP1,J)
P4I=P4I+P4
AIJ(I,J,IJ_UJET)=AIJ(I,J,IJ_UJET)+U(I,J,JET)
AIJ(I,J,IJ_VJET)=AIJ(I,J,IJ_VJET)+V(I,J,JET)
I=IP1
END DO
APJ(J,2)=APJ(J,2)+P4I
DO L=1,LM
PU4I=0.
PV4I=0.
PUV4I=0.
I=IM
DO IP1=1,IM
P4=PLIJ(L,I,J-1)+PLIJ(L,IP1,J-1)+PLIJ(L,I,J)+PLIJ(L,IP1,J)
IF(L.EQ.LS1) P4I=FIM*P4
PU4I=PU4I+P4*U(I,J,L)
PV4I=PV4I+P4*V(I,J,L)
PUV4I=PUV4I+P4*U(I,J,L)*V(I,J,L)
T4=TX(I,J-1,L)+TX(IP1,J-1,L)+TX(I,J,L)+TX(IP1,J,L)
Z4=PHI(I,J-1,L)+PHI(IP1,J-1,L)+PHI(I,J,L)+PHI(IP1,J,L)
AIJ(I,J,IJ_DSEV)=AIJ(I,J,IJ_DSEV)+P4*(SHA*T4+Z4)*V(I,J,L)
* *DSIG(L)*DXV(J)
SP2=PLIJ(L,IP1,J-1)+PLIJ(L,IP1,J)
AIJ(IP1,J,IJ_PVQ)=AIJ(IP1,J,IJ_PVQ)+SP2
* *( I=IP1
END DO
AJL(J,L,JL_ZMFNTMOM)=AJL(J,L,JL_ZMFNTMOM)+PU4I*PV4I/P4I
AJL(J,L,JL_TOTNTMOM)=AJL(J,L,JL_TOTNTMOM)+PUV4I
END DO
END DO
C****
C**** EVEN LEVEL GEOPOTENTIALS, VERTICAL WINDS AND VERTICAL TRANSPORTS
C****
DO L=1,LM-1
DO J=J_0,J_1
DO I=1,IMAXJ(J)
PIJ=PLIJ(L,I,J)
PE=SIGE(L+1)*PIJ+PTOP
PKE=PE**KAPA
THETA=THBAR( W(I,J,L)=SD(I,J,L)*THETA*PKE/PE
END DO
END DO
END DO
C****
C**** AVAILABLE POTENTIAL ENERGY
C****
C**** SET UP FOR CALCULATION
DO 710 L=1,LM
710 GMEAN(L)=0.
DO 740 J=J_0,J_1
DO 720 I=1,IMAXJ(J)
720 SQRTP(I)=SQRT(C**** GMEAN CALCULATED FOR EACH LAYER, THJL, THSQJL ARRAYS FILLED
DO 730 L=1,LM
LDN=LDNA(L)
LUP=LUPA(L)
THJL(J,L)=0.
THSQJL(J,L)=0.
DO 730 I=1,IMAXJ(J)
THJL(J,L)=THJL(J,L)+T(I,J,L)*SQRTP(I)
THSQJL(J,L)=THSQJL(J,L)+T(I,J,L)*T(I,J,L)*P(I,J)
730 GMEAN(L)=GMEAN(L)+(SIG(L)*P(I,J)+PTOP)*( * DXYP(J)/( 740 CONTINUE
C**** CALCULATE APE
DO 760 L=1,LM
LP1=LUPA(L)
LM1=LDNA(L)
IF(GRID%HAVE_SOUTH_POLE) THEN
THJL(1,L)=THJL(1,L)*FIM
THSQJL(1,L)=THSQJL(1,L)*FIM
ENDIF
IF(GRID%HAVE_NORTH_POLE) THEN
THJL(JM,L)=THJL(JM,L)*FIM
THSQJL(JM,L)=THSQJL(JM,L)*FIM
ENDIF
THGM=0.
DO 750 J=J_0,J_1
750 THGM=THGM+THJL(J,L)*DXYP(J)
THGM=THGM/AREAG
GMEANL=GMEAN(L)/((SIG(LM1)-SIG(LP1))*AREAG)
DO 760 J=J_0,J_1
760 AJL(J,L,JL_APE)=AJL(J,L,JL_APE)+
& (THSQJL(J,L)-2.*THJL(J,L)*THGM+THGM*THGM*
* FIM)/GMEANL
C****
C**** CERTAIN HORIZONTAL WIND AVERAGES
C****
DO L=1,LM
DO J=J_0STG,J_1STG
DO I=I135W,I110W ! EAST PACIFIC
AJL(J,L,JL_UEPAC)=AJL(J,L,JL_UEPAC)+U(I,J,L)
AJL(J,L,JL_VEPAC)=AJL(J,L,JL_VEPAC)+V(I,J,L)
END DO
DO I=I150E,IM ! WEST PACIFIC
AJL(J,L,JL_UWPAC)=AJL(J,L,JL_UWPAC)+U(I,J,L)
AJL(J,L,JL_VWPAC)=AJL(J,L,JL_VWPAC)+V(I,J,L)
END DO
END DO
DO J=MAX(J_0,J5SUV),MIN(J_1,J5NUV)
DO I=1,IM
AIL(I,L,IL_UEQ)=AIL(I,L,IL_UEQ)+U(I,J,L)
AIL(I,L,IL_VEQ)=AIL(I,L,IL_VEQ)+V(I,J,L)
END DO
END DO
DO J=MAX(J_0,J5S),MIN(J_1,J5N)
DO I=1,IM
AIL(I,L,IL_TEQ)=AIL(I,L,IL_TEQ)+(TX(I,J,L)-TF)
AIL(I,L,IL_QEQ)=AIL(I,L,IL_QEQ)+Q(I,J,L)/QSAT(TX(I,J,L),LHE
* ,PMID(L,I,J))
END DO
END DO
IF(J_0 <= J50N .and. J_1 >= J50N) THEN
DO I=1,IM
AIL(I,L,IL_T50N)=AIL(I,L,IL_T50N)+(TX(I,J50N,L)-TF)
AIL(I,L,IL_U50N)=AIL(I,L,IL_U50N)+( * U(I,J50N+1,L))
END DO
ENDIF
IF(J_0 <= J70N .and. J_1 >= J70N) THEN
DO I=1,IM
AIL(I,L,IL_T70N)=AIL(I,L,IL_T70N)+(TX(I,J70N,L)-TF)
AIL(I,L,IL_U70N)=AIL(I,L,IL_U70N)+( * U(I,J70N+1,L))
END DO
ENDIF
END DO
C****
C**** CERTAIN VERTICAL WIND AVERAGES
C****
DO L=1,LM-1
DO J=J_0S,J_1S
DO I=I135W,I110W ! EAST PACIFIC
AJL(J,L,JL_WEPAC)=AJL(J,L,JL_WEPAC)+W(I,J,L)
END DO
DO I=I150E,IM ! WEST PACIFIC
AJL(J,L,JL_WWPAC)=AJL(J,L,JL_WWPAC)+W(I,J,L)
END DO
END DO
DO I=1,IM
DO J=MAX(J_0,J5S),MIN(J_1,J5N) ! +/- 5 DEG (APPROX.)
AIL(I,L,IL_WEQ) =AIL(I,L,IL_WEQ)+W(I,J,L)
END DO
IF(J_0 <= J50N .and. J_1 >= J50N)
& AIL(I,L,IL_W50N)=AIL(I,L,IL_W50N)+W(I,J50N,L)
IF(J_0 <= J70N .and. J_1 >= J70N)
& AIL(I,L,IL_W70N)=AIL(I,L,IL_W70N)+W(I,J70N,L)
END DO
END DO
C****
C**** ELIASSEN PALM FLUX
C****
C**** NORTHWARD TRANSPORT
!Not necessary here, done above CALL CHECKSUM(grid, P, 752, "DIAG.f")
!Not necessary here, done above CALL HALO_UPDATE(grid, P, FROM=SOUTH)
CALL CHECKSUM(grid, DXYN, 754, "DIAG.f")
CALL HALO_UPDATE(grid, DXYN, FROM=SOUTH)
CALL CHECKSUM(grid, T, 756, "DIAG.f")
CALL HALO_UPDATE(grid, T, FROM=SOUTH)
DO 868 J=J_0STG,J_1STG
I=IM
DO 862 IP1=1,IM
PDA(I)=.5*((P(I,J)+P(IP1,J))*DXYS(J)+( * DXYN(J-1))
PSEC(I)=PDA(I)*BYDXYV(J)
862 I=IP1
DO 868 L=1,LM
DUDP=0.
DTHDP=0.
UMN=0.
THMN=0.
LDN=LDNA(L)
LUP=LUPA(L)
I=IM
DO 864 IP1=1,IM
DUDP=DUDP+U(I,J,LUP)-U(I,J,LDN)
DTHDP=DTHDP+T(I,J,LUP)+T(I,J-1,LUP)-T(I,J,LDN)-T(I,J-1,LDN)
UMN=UMN+U(I,J,L)
THMN=THMN+T(I,J,L)+T(I,J-1,L)
THSEC(I)=T(I,J,L)+T(IP1,J,L)+T(I,J-1,L)+T(IP1,J-1,L)
864 I=IP1
UMN=UMN*BYIM
THMN=2.*THMN/FIM
FPHI=0.
SMALL=.0002*FIM*T(1,J,L)
c IF (DTHDP.LT.SMALL) WRITE (6,999) J,L,DTHDP,SMALL
IF (DTHDP.LT.SMALL) DTHDP=SMALL
DO 866 I=1,IM
SP=PSEC(I)
IF(L.GE.LS1) SP=PSFMPT
866 FPHI=FPHI+SP*V(I,J,L)*(.5*(THSEC(I)-THMN)*DUDP/DTHDP
* -U(I,J,L)+UMN)
868 AJL(J,L,JL_EPFLXN)=AJL(J,L,JL_EPFLXN)+FPHI
C**** VERTICAL TRANSPORT
!Not necessary here, done above CALL CHECKSUM(grid, U, 794, "DIAG.f")
!Not necessary here, done above CALL HALO_UPDATE(grid, U, FROM=NORTH)
CALL CHECKSUM(grid, V, 796, "DIAG.f")
CALL HALO_UPDATE(grid, V, FROM=NORTH)
DO 878 J=J_0S,J_1S
PITMN=0.
DO 870 I=1,IM
870 PITMN=PITMN+PIT(I,J)
PITMN=PITMN/FIM
DO 878 L=1,LM-1
IF(L.GE.LS1-1) PITMN=0.
THMN=0.
SDMN=0.
DTHDP=0.
DO 872 I=1,IM
DTHDP=DTHDP+T(I,J,L+1)-T(I,J,L)
THMN=THMN+T(I,J,L+1)+T(I,J,L)
872 SDMN=SDMN+SD(I,J,L)
SMALL=.0001*FIM*T(1,J,L+1)
c IF (DTHDP.LT.SMALL) WRITE (6,999) J,L,DTHDP,SMALL
IF (DTHDP.LT.SMALL) DTHDP=SMALL
THMN=THMN/FIM
SDMN=SDMN/FIM
DUDX=0.
PVTHP=0.
SDPU=0.
IM1=IM
DO 874 I=1,IM
DUDX=DUDX+DXV(J+1)*( * ( UPE=U(IM1,J,L)+U(IM1,J+1,L)+U(I,J,L)+U(I,J+1,L)+
* U(IM1,J,L+1)+U(IM1,J+1,L+1)+U(I,J,L+1)+U(I,J+1,L+1)
VPE=V(IM1,J,L)+V(IM1,J+1,L)+V(I,J,L)+V(I,J+1,L)+
* V(IM1,J,L+1)+V(IM1,J+1,L+1)+V(I,J,L+1)+V(I,J+1,L+1)
DP=(SIG(L)-SIG(L+1))*P(I,J)
IF(L.GE.LS1) DP=(SIG(L)-SIG(L+1))*PSFMPT
IF(L.EQ.LS1-1) DP=P(I,J)*SIG(L)-PSFMPT*SIG(LS1)
PVTHP=PVTHP+DP*VPE*( PITIJ=PIT(I,J)
IF(L.GE.LS1-1) PITIJ=0.
SDPU=SDPU+(SD(I,J,L)-SDMN+(PITIJ-PITMN)*SIGE(L+1))*UPE
874 IM1=I
AJL(J,L,JL_EPFLXV)=AJL(J,L,JL_EPFLXV)+
& (.5*FIM*FCOR(J)-.25*DUDX)*PVTHP/DTHDP + SDPU
878 CONTINUE
C**** ACCUMULATE TIME USED IN DIAGA
CALL TIMEOUT(MBEGIN,MDIAG,MDYN)
RETURN
999 FORMAT (' DTHETA/DP IS TOO SMALL AT J=',I4,' L=',I4,2F15.6)
ENTRY DIAGA0 1
C****
C**** INITIALIZE TJL0 ARRAY (FROM PRIOR TO ADVECTION)
C****
CALL GET(grid, J_STRT=J_0, J_STOP=J_1)
DO L=1,LM
DO J=J_0,J_1
THI=0.
DO I=1,IMAXJ(J)
THI=THI+T(I,J,L)
END DO
TJL0(J,L)=THI
END DO
END DO
RETURN
C****
END SUBROUTINE DIAGA
SUBROUTINE DIAGB 1,18
!@sum DIAGB calculate constant pressure diagnostics from within DYNAM
C****
C**** CONTENTS OF AJK(J,K,N) (SUM OVER LONGITUDE AND TIME OF)
C**** See jks_defs for contents
C****
C**** CONTENTS OF AIJK(I,J,K,N) (SUM OVER TIME OF)
C**** See ijks_defs for contents
C****
USE CONSTANT, only : lhe,omega,sha,tf,teeny
USE MODEL_COM, only :
& im,imh,fim,byim,jm,jeq,lm,ls1,idacc,ptop,psfmpt,jdate,
& mdyn,mdiag, ndaa,sig,sige,dsig,Jhour,u,v,t,p,q,wm,km=>lm
USE GEOM, only :
& COSV,DXV,DXYN,DXYP,DXYS,DXYV,DYP,DYV,FCOR,IMAXJ,RADIUS
USE DAGCOM, only : ajk,aijk,speca,adiurn,nspher,hdiurn,
& nwav_dag,ndiupt,hr_in_day,ijk_u,ijk_v,ijk_t,ijk_q,ijk_dp
* ,ijk_dse,klayer,idd_w,ijdd,
& JK_DPA,JK_DPB,JK_TEMP,JK_HGHT,JK_Q,JK_THETA,
& JK_RH,JK_U,JK_V,JK_ZMFKE,JK_TOTKE,JK_ZMFNTSH,
& JK_TOTNTSH,JK_ZMFNTGEO,JK_TOTNTGEO,JK_ZMFNTLH,
& JK_TOTNTLH,JK_ZMFNTKE,JK_TOTNTKE,JK_ZMFNTMOM,JK_TOTNTMOM,
& JK_P2KEDPGF,JK_DPSQR,JK_NPTSAVG,
& JK_VVEL,JK_ZMFVTDSE,JK_TOTVTDSE,JK_ZMFVTLH,JK_TOTVTLH,
& JK_VTGEOEDDY,JK_BAREKEGEN,JK_POTVORT,JK_VTPV,
& JK_VTPVEDDY,JK_NPTSAVG1,JK_TOTVTKE,
& JK_VTAMEDDY,JK_TOTVTAM,JK_SHETH,JK_DUDTMADV,JK_DTDTMADV,
& JK_DUDTTEM,JK_DTDTTEM,JK_EPFLXNCP,JK_EPFLXVCP,
& JK_UINST,JK_TOTDUDT,JK_TINST,
& JK_TOTDTDT,JK_EDDVTPT,JK_CLDH2O
USE DYNAMICS, only : phi,dut,dvt,plij
USE DIAG_LOC, only : w,tx,pm,pl,pmo,plo
USE DOMAIN_DECOMP, only : GET, CHECKSUM, HALO_UPDATE, GRID
USE DOMAIN_DECOMP, only : SOUTH, NORTH
IMPLICIT NONE
REAL*8, DIMENSION(IMH+1,NSPHER) :: KE
REAL*8, DIMENSION(IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM) ::
& ZX,STB
REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM) ::
& STJK,DPJK,UJK,VJK,WJK,TJK,
& PSIJK,UP,TY,PSIP,WTJK,UVJK,WUJK
REAL*8, DIMENSION(IM) :: PSEC,X1
REAL*8, DIMENSION(LM) :: SHETH,DPM,DTH
INTEGER ::
& I,IH,IHM,IM1,INCH,INCHM,IP1,IZERO,J,J45N,
& JHEMI,K,KDN,KR,KS,KS1,KSPHER,KUP,KX,L,
& LUP,MBEGIN,N,NM
REAL*8 ::
& BYDP,BYFIM,DP,DPDN,DPDX,DPDY,
& DPE,DPI,DPK,DPSQI,DPUP,DPUV,DUTI,DUTK,DVTI,DVTK,FIMI,
& PAI,PAK,PDN,PHIPI,PMK,PQ4I,PQ4K,PQV4I,PS,PS4I,
& PS4K,PSIY,PSV4I,PT4I,PT4K,PTK,PTV4I,PUI,PUK,PUP,
& PUVI,PV2,PV2I,PVI,PVK,PWWI,PWWVI,PY,PZ4I,PZ4K,
& PZV4I,QK,QKI,QLH,QPI,QSATL,RHPI,
& SMALL,SP,SP2,SQRTDP,THK,THKI,THPI,TK,TKI,TPI,
& UDUTI,UDX,UEARTH,UK,UKI,UY,VDVTI,VK,VSTAR,W2,W2I,W4,
& W4I,WI,WKE4I,WMPI,WNP,WPA2I,WPV4I,WQI,WSP,WSTAR,WTHI,
& WTI,WU4I,WUP,WZI,ZK,ZKI
REAL*8, PARAMETER :: BIG=1.E20
REAL*8 :: QSAT
INTEGER :: J_0, J_1, J_0S, J_1S, J_0STG, J_1STG
CALL GETTIME(MBEGIN)
CALL GET(grid, J_STRT=J_0, J_STOP=J_1,
& J_STRT_SKP=J_0S, J_STOP_SKP=J_1S,
& J_STRT_STGR=J_0STG, J_STOP_STGR=J_1STG)
C****
C**** INTERNAL QUANTITIES T,TH,Q,RH
C****
QLH=LHE
DO 170 J=J_0,J_1
DO 170 K=1,KM
DPI=0.
TPI=0.
PHIPI=0.
QPI=0.
WMPI=0.
THPI=0.
RHPI=0.
FIMI=0.
DO 160 I=1,IMAXJ(J)
C**** FIND L=L(K) AND LUP=L(K+1) S.T. P(LUP).GT.P(K+1)
SP=PLIJ(K,I,J)
PS=SP+PTOP
IF (PM(K+1).GE.PS) GO TO 160
L=1
PDN=PS
IF (PM(K).GE.PS) GO TO 120
PDN=PM(K)
110 IF (PM(K).GT.SP*SIGE(L+1)+PTOP) GO TO 120
L=L+1
GO TO 110
120 LUP=L
130 IF (PM(K+1).GE.SP*SIGE(LUP+1)+PTOP) GO TO 140
LUP=LUP+1
GO TO 130
140 CONTINUE
C**** ACCUMULATE HERE
DPI=DPI+PDN-PM(K+1)
FIMI=FIMI+1.
150 PUP=SP*SIGE(L+1)+PTOP
IF (LUP.EQ.L) PUP=PM(K+1)
DP=PDN-PUP
TPI=TPI+(TX(I,J,L)-TF)*DP
PHIPI=PHIPI+PHI(I,J,L)*DP
QPI=QPI+Q(I,J,L)*DP
WMPI=WMPI+WM(I,J,L)*DP
CW IF(WMPI.GT.1.E-3) WRITE(6,169) I,J,L,DP,WM(I,J,L),WMPI
CW169 FORMAT(1X,'1616--',3I5,3E15.2)
THPI=THPI+T(I,J,L)*DP
QSATL=QSAT(TX(I,J,L),QLH,SIG(L)*SP+PTOP)
IF (QSATL.GT.1.) QSATL=1.
RHPI=RHPI+Q(I,J,L)*DP/QSATL
IF (L.EQ.LUP) GO TO 160
L=L+1
PDN=SP*SIGE(L)+PTOP
GO TO 150
160 CONTINUE
AJK(J,K,JK_NPTSAVG1)=AJK(J,K,JK_NPTSAVG1)+FIMI
AJK(J,K,JK_DPA)=AJK(J,K,JK_DPA)+DPI
AJK(J,K,JK_TEMP)=AJK(J,K,JK_TEMP)+TPI
AJK(J,K,JK_HGHT)=AJK(J,K,JK_HGHT)+PHIPI
AJK(J,K,JK_Q)=AJK(J,K,JK_Q)+QPI
AJK(J,K,JK_THETA)=AJK(J,K,JK_THETA)+THPI
AJK(J,K,JK_RH)=AJK(J,K,JK_RH)+RHPI
AJK(J,K,JK_CLDH2O)=AJK(J,K,JK_CLDH2O)+WMPI
TJK(J,K)=THPI/(DPI+teeny)
IF (IDACC(4).EQ.1) AJK(J,K,JK_TINST)=TJK(J,K)
AJK(J,K,JK_TOTDTDT)=TJK(J,K)-AJK(J,K,JK_TINST)
170 CONTINUE
C****
C**** CALCULATE STABILITY AT ODD LEVELS ON PU GRID
C****
DO 230 J=J_0,J_1
I=IMAXJ(J)
DO 230 IP1=1,IMAXJ(J)
SP2=P(I,J)+P(IP1,J)
SP=.5*SP2
DO 175 L=1,LS1-1
PLO(L)=SP*SIG(L)+PTOP
175 PL(L)=SP*SIGE(L)+PTOP
DO 180 L=1,LM-1
DTH(L)=( * (2.*(PLO(L)-PLO(L+1)))
180 CONTINUE
DO 220 K=1,KM
STB(I,J,K)=0.
IF (PM(K+1).GE.PL(1)) GO TO 220
PMK=PMO(K)
IF (PM(K).GT.PL(1)) PMK=.5*(SP+PTOP+PM(K+1))
L=2
IF (PMK.GE.PL(2)) GO TO 210
190 LUP=L+1
IF (L.EQ.LM) GO TO 210
IF (PMK.GE.PL(LUP)) GO TO 200
L=LUP
GO TO 190
200 DPUP=PMK-PL(LUP)
DPDN=PL(L)-PMK
STB(I,J,K)=(DTH(L-1)*DPUP+DTH(L)*DPDN)/(DPUP+DPDN+teeny)
GO TO 220
C**** SPECIAL CASES, L=2, L=LM
210 STB(I,J,K)=DTH(L-1)
220 CONTINUE
230 I=IP1
C**** CALCULATE STJK; THE MEAN STATIC STABILITY
DO 260 J=J_0,J_1
DO 260 K=1,KM
STJK(J,K)=0.
DPJK(J,K)=0.
I=IMAXJ(J)
DO 250 IP1=1,IMAXJ(J)
PS=.5*( IF (PM(K+1).GT.PS) GO TO 250
STJK(J,K)=STJK(J,K)+STB(I,J,K)
DPJK(J,K)=DPJK(J,K)+1.
250 I=IP1
STJK(J,K)=STJK(J,K)/(DPJK(J,K)+teeny)
SMALL=.0001
IF (ABS(STJK(J,K)).LT.SMALL) STJK(J,K)=-SMALL
260 CONTINUE
C****
C**** CONSTANT PRESSURE DIAGNOSTICS: FLUX, ENERGY, ANGULAR MOMENTUM
C****
IF(GRID%HAVE_SOUTH_POLE) THEN
ZX(:,1,:)=0.
ENDIF
C Needs to check later if all these halo calls are necessary.
C DIAGA may have contained relevant halo calls
C and since DIAGB is called immediately after DIAGA
C there may not be a need for these calls if
C the concerned arrays have not been updated
C from the previous halo call.
CALL CHECKSUM(grid, P, 1064, "DIAG.f")
CALL HALO_UPDATE(grid, P, FROM=SOUTH)
CALL CHECKSUM(grid, TX, 1066, "DIAG.f")
CALL HALO_UPDATE(grid, TX, FROM=SOUTH)
CALL CHECKSUM(grid, PHI, 1068, "DIAG.f")
CALL HALO_UPDATE(grid, PHI, FROM=SOUTH)
CALL CHECKSUM(grid, Q, 1070, "DIAG.f")
CALL HALO_UPDATE(grid, Q, FROM=SOUTH)
CALL CHECKSUM(grid, T, 1072, "DIAG.f")
CALL HALO_UPDATE(grid, T, FROM=SOUTH)
CALL CHECKSUM(grid, ZX, 1074, "DIAG.f")
CALL HALO_UPDATE(grid, ZX, FROM=SOUTH)
CALL CHECKSUM(grid, STJK, 1076, "DIAG.f")
CALL HALO_UPDATE(grid, STJK, FROM=SOUTH)
CALL CHECKSUM(grid, DXYN, 1078, "DIAG.f")
CALL HALO_UPDATE(grid, DXYN, FROM=SOUTH)
DO 390 J=J_0STG,J_1STG
I=IM
DO 280 IP1=1,IM
PSEC(I)=.25*( DO 270 K=1,KM
270 ZX(I,J,K)=0.
DO 275 L=1,LS1-1
DUT(I,J,L)=DUT(I,J,L)/(PSEC(I)*DXYV(J)*DSIG(L))
275 DVT(I,J,L)=DVT(I,J,L)/(PSEC(I)*DXYV(J)*DSIG(L))
DO 276 L=LS1,LM
DUT(I,J,L)=DUT(I,J,L)/(PSFMPT*DXYV(J)*DSIG(L))
276 DVT(I,J,L)=DVT(I,J,L)/(PSFMPT*DXYV(J)*DSIG(L))
280 I=IP1
DO 350 K=1,KM
DPI=0.
DPSQI=0.
FIMI=0.
PUI=0.
PVI=0.
PWWI=0.
PT4I=0.
PTV4I=0.
PZ4I=0.
PZV4I=0.
PQ4I=0.
PQV4I=0.
PWWVI=0.
PUVI=0.
DVTI=0.
VDVTI=0.
DUTI=0.
UDUTI=0.
PS4I=0.
PSV4I=0.
I=IM
DO 340 IP1=1,IM
SP=PSEC(I)
PS=SP+PTOP
DO 286 L=1,LS1-1
286 PL(L)=SP*SIGE(L)+PTOP
IF (PM(K+1).GE.PS) GO TO 336
L=1
PDN=PS
IF (PM(K).GE.PS) GO TO 300
PDN=PM(K)
290 IF (PM(K).GT.PL(L+1)) GO TO 300
L=L+1
GO TO 290
300 LUP=L
310 IF (PM(K+1).GE.PL(LUP+1)) GO TO 320
LUP=LUP+1
GO TO 310
320 CONTINUE
DPK=PDN-PM(K+1)
PUK=0.
PVK=0.
PT4K=0.
PZ4K=0.
PQ4K=0.
DUTK=0.
DVTK=0.
PS4K=0.
C**** INTERPOLATE HERE
330 PUP=PL(L+1)
IF (LUP.EQ.L) PUP=PM(K+1)
DP=PDN-PUP
PUK=PUK+DP*U(I,J,L)
PVK=PVK+DP*V(I,J,L)
PT4K=PT4K+DP*(TX(I,J-1,L)+TX(IP1,J-1,L)+TX(I,J,L)+TX(IP1,J,L))
PZ4K=PZ4K+DP*(PHI(I,J-1,L)+PHI(IP1,J-1,L)+PHI(I,J,L)+PHI(IP1,J,L))
PQ4K=PQ4K+DP*( DUTK=DUTK+DP*DUT(I,J,L)
DVTK=DVTK+DP*DVT(I,J,L)
PS4K=PS4K+DP*( IF (LUP.EQ.L) GO TO 332
L=L+1
PDN=PL(L)
GO TO 330
C**** ACCUMULATE HERE
332 FIMI=FIMI+1.
DPI=DPI+DPK
DPSQI=DPSQI+DPK*DPK
IF (DPK.LT.teeny) DPK=teeny
BYDP=1./DPK
PUI=PUI+PUK
PVI=PVI+PVK
PWWI=PWWI+BYDP*(PUK*PUK+PVK*PVK)
PWWVI=PWWVI+BYDP*BYDP*(PUK*PUK+PVK*PVK)*PVK
PUVI=PUVI+BYDP*PUK*PVK
PT4I=PT4I+PT4K
PTV4I=PTV4I+BYDP*PT4K*PVK
PZ4I=PZ4I+PZ4K
PZV4I=PZV4I+BYDP*PZ4K*PVK
PQ4I=PQ4I+PQ4K
PQV4I=PQV4I+BYDP*PQ4K*PVK
DVTI=DVTI+DVTK
VDVTI=VDVTI+BYDP*PVK*DVTK
DUTI=DUTI+DUTK
UDUTI=UDUTI+BYDP*PUK*DUTK
!! IF(SKIPSE.EQ.1.) GO TO 334
AIJK(I,J,K,IJK_U) =AIJK(I,J,K,IJK_U) +PUK
AIJK(I,J,K,IJK_V) =AIJK(I,J,K,IJK_V) +PVK
AIJK(I,J,K,IJK_DSE)=AIJK(I,J,K,IJK_DSE)+SHA*PT4K+PZ4K
AIJK(I,J,K,IJK_DP) =AIJK(I,J,K,IJK_DP) +DPK
AIJK(I,J,K,IJK_T) =AIJK(I,J,K,IJK_T) +PT4K
AIJK(I,J,K,IJK_Q) =AIJK(I,J,K,IJK_Q) +PQ4K
C**** EDDY TRANSPORT OF THETA; VORTICITY
334 PS4I=PS4I+PS4K
PSV4I=PSV4I+BYDP*PVK*PS4K
UDX=BYDP*PUK*DXV(J)
ZX(I,J,K)=-UDX
IF (ZX(I,J-1,K).LT.BIG) ZX(I,J-1,K)=ZX(I,J-1,K)+UDX
IF (ZX(I,J-1,K).GE.BIG) ZX(I,J-1,K)=0.
GO TO 340
336 ZX(I,J,K)=BIG
ZX(I,J-1,K)=0.
340 I=IP1
DPM(K)=DPI/(FIMI+teeny)
DPJK(J,K)=DPI
AJK(J,K,JK_DPB)=AJK(J,K,JK_DPB)+DPI
AJK(J,K,JK_DPSQR)=AJK(J,K,JK_DPSQR)+DPSQI
AJK(J,K,JK_NPTSAVG)=AJK(J,K,JK_NPTSAVG)+FIMI
IF (DPI.LT.teeny) DPI=teeny
AJK(J,K,JK_U)=AJK(J,K,JK_U)+PUI
AJK(J,K,JK_V)=AJK(J,K,JK_V)+PVI
AJK(J,K,JK_ZMFKE)=AJK(J,K,JK_ZMFKE)+(PUI*PUI+PVI*PVI)/DPI
AJK(J,K,JK_TOTKE)=AJK(J,K,JK_TOTKE)+PWWI
AJK(J,K,JK_ZMFNTSH)=AJK(J,K,JK_ZMFNTSH)+PT4I*PVI/DPI
AJK(J,K,JK_TOTNTSH)=AJK(J,K,JK_TOTNTSH)+PTV4I
AJK(J,K,JK_ZMFNTGEO)=AJK(J,K,JK_ZMFNTGEO)+PZ4I*PVI/DPI
AJK(J,K,JK_TOTNTGEO)=AJK(J,K,JK_TOTNTGEO)+PZV4I
AJK(J,K,JK_ZMFNTLH)=AJK(J,K,JK_ZMFNTLH)+PQ4I*PVI/DPI
AJK(J,K,JK_TOTNTLH)=AJK(J,K,JK_TOTNTLH)+PQV4I
AJK(J,K,JK_ZMFNTKE)=AJK(J,K,JK_ZMFNTKE)+PWWI*PVI/DPI
AJK(J,K,JK_TOTNTKE)=AJK(J,K,JK_TOTNTKE)+PWWVI
AJK(J,K,JK_ZMFNTMOM)=AJK(J,K,JK_ZMFNTMOM)+PUI*PVI/DPI
AJK(J,K,JK_TOTNTMOM)=AJK(J,K,JK_TOTNTMOM)+PUVI
AJK(J,K,JK_P2KEDPGF)=AJK(J,K,JK_P2KEDPGF)+VDVTI+UDUTI-
* (PUI*DUTI+PVI*DVTI)/DPI
SHETH(K)=(PSV4I-PS4I*PVI/DPI)*DXYV(J)/(STJK(J-1,K)*DXYN(J-1)+
* STJK(J,K)*DXYS(J))
UJK(J,K)=PUI/DPI
VJK(J,K)=PVI/DPI
PSIJK(J,K)=.25*SHETH(K)/DPI
UVJK(J,K)=(PUVI-PUI*PVI/DPI)/DPI
IF (IDACC(4).EQ.1) AJK(J,K,JK_UINST)=UJK(J,K)
AJK(J,K,JK_TOTDUDT)=UJK(J,K)-AJK(J,K,JK_UINST)
350 AJK(J,K,JK_SHETH)=AJK(J,K,JK_SHETH)+SHETH(K)
390 CONTINUE
C****
C**** VERTICAL MASS FLUXES W(I,J,K)
C****
DO 400 I=1,IM
DO 400 J=J_0,J_1
DO 400 K=1,KM
DUT(I,J,K)=0.
DVT(I,J,K)=0.
400 W(I,J,K)=0.
CALL CHECKSUM(grid, U, 1240, "DIAG.f")
CALL HALO_UPDATE(grid, U, FROM=NORTH)
C**** EASTWARD MASS FLUX DUT (PU POINTS)
DO 460 J=J_0S,J_1S
DO 460 K=1,KM
I=IM
DO 460 IP1=1,IM
SP=.5*( DO 405 L=1,LS1-1
405 PL(L)=SP*SIGE(L)+PTOP
IF (PM(K+1).GE.SP+PTOP) GO TO 460
L=1
PDN=SP+PTOP
IF (PM(K).GE.SP+PTOP) GO TO 420
PDN=PM(K)
410 IF (PM(K).GT.PL(L+1)) GO TO 420
L=L+1
GO TO 410
420 LUP=L
430 IF (PM(K+1).GE.PL(LUP+1)) GO TO 440
LUP=LUP+1
GO TO 430
440 CONTINUE
C**** CALCULATE HERE
450 PUP=PL(L+1)
IF (LUP.EQ.L) PUP=PM(K+1)
DPDY=(PDN-PUP)*DYP(3)
DUT(I,J,K)=DUT(I,J,K)+DPDY*( IF (LUP.EQ.L) GO TO 460
L=L+1
PDN=PL(L)
GO TO 450
460 I=IP1
C**** NORTHWARD MASS FLUX DVT (PV POINTS)
C P already halo'ed; no need CALL CHECKSUM(grid, P, __LINE__, __FILE__)
C P already halo'ed; no need CALL HALO_UPDATE(grid, P, FROM=SOUTH)
DO 520 J=J_0STG,J_1STG
DO 520 K=1,KM
IM1=IM
DO 520 I=1,IM
SP=.5*( DO 465 L=1,LS1-1
465 PL(L)=SP*SIGE(L)+PTOP
IF (PM(K+1).GE.SP+PTOP) GO TO 520
L=1
PDN=SP+PTOP
IF (PM(K).GE.SP+PTOP) GO TO 480
PDN=PM(K)
470 IF (PM(K).GT.PL(L+1)) GO TO 480
L=L+1
GO TO 470
480 LUP=L
490 IF (PM(K+1).GE.PL(LUP+1)) GO TO 500
LUP=LUP+1
GO TO 490
500 CONTINUE
C**** CALCULATE HERE
510 PUP=PL(L+1)
IF (LUP.EQ.L) PUP=PM(K+1)
DPDX=(PDN-PUP)*DXV(J)
DVT(I,J,K)=DVT(I,J,K)+DPDX*( IF (LUP.EQ.L) GO TO 520
L=L+1
PDN=PL(L)
GO TO 510
520 IM1=I
CALL CHECKSUM(grid, DVT, 1309, "DIAG.f")
CALL HALO_UPDATE(grid, DVT, FROM=NORTH)
DO 560 K=KM,1,-1
C**** POLAR VERTICAL MASS FLUX
IF(GRID%HAVE_SOUTH_POLE) THEN
W(1,1,K)=0.
IF (K.LT.KM) W(1,1,K)=W(1,1,K+1)
DO I=1,IM
W(1,1,K)=W(1,1,K)-.5*DVT(I,2,K)
ENDDO
ENDIF
IF(GRID%HAVE_NORTH_POLE) THEN
W(1,JM,K)=0.
IF (K.LT.KM) W(1,JM,K)=W(1,JM,K+1)
DO I=1,IM
W(1,JM,K)=W(1,JM,K)+.5*DVT(I,JM,K)
ENDDO
ENDIF
C**** NON-POLAR VERTICAL MASS FLUX
WUP=0.
DO 560 J=J_0S,J_1S
IM1=IM
DO 560 I=1,IM
IF (K.LT.KM) WUP=W(I,J,K+1)
W(I,J,K)=WUP+.5*(DUT(IM1,J,K)-DUT(I,J,K)+
* DVT(I,J,K)-DVT(I,J+1,K))
560 IM1=I
C**** ZERO OUT SUBSURFACE VERTICAL WINDS
DO J=J_0,J_1
DO I=1,IM
PS=P(I,J)+PTOP
K=2
DO WHILE(PM(K+1).GE.PS)
W(I,J,K)=0.
K=K+1
ENDDO
ENDDO
ENDDO
C**** ACCUMULATE ALL VERTICAL WINDS
DO 558 J=J_0,J_1
DO 558 I=1,IM
DO KR=1,NDIUPT
IF(I.EQ.IJDD(1,KR).AND.J.EQ.IJDD(2,KR)) THEN
C**** Warning: This diagnostic has 3 flaws (?)
C**** 1 - It assumes that DTsrc=1hr, (DTsrc=3600.)
C**** 2 - since DTdaa-Ndaa*DTsrc=2*DTdyn rather than 0,
C**** some hours are skipped once in a while
C**** 3 - Some of the first Ndaa hours are skipped at the
C**** beginning of a month and overcounted at the end;
C**** this happens to balance out, if and only if
C**** mod(days_in_month,ndaa)=0 (i.e. February if Ndaa=7)
IH=JHOUR+1
IHM = IH+(JDATE-1)*24
DO INCH=1,NDAA
IF(IH.GT.HR_IN_DAY) IH=IH-HR_IN_DAY
ADIURN(IH,IDD_W,KR)=ADIURN(IH,IDD_W,KR)+1.E5*W(I,J,3)
* /DXYP(J)
HDIURN(IHM,IDD_W,KR)=HDIURN(IHM,IDD_W,KR)+1.E5*W(I,J,3)
* /DXYP(J)
IH=IH+1
IHM=IHM+1
END DO
END IF
END DO
558 CONTINUE
DO 565 J=J_0,J_1
DO 565 K=1,KM
WI=0.
DO 562 I=1,IMAXJ(J)
562 WI=WI+W(I,J,K)
565 AJK(J,K,JK_VVEL)=AJK(J,K,JK_VVEL)+WI
C****
C**** ACCUMULATE T,Z,Q VERTICAL TRANSPORTS
C****
DO 610 J=J_0,J_1
DO 610 K=2,KM
WI=0.
TKI=0.
QKI=0.
ZKI=0.
WTI=0.
WQI=0.
WZI=0.
THKI=0.
WTHI=0.
FIMI=0.
DO 600 I=1,IMAXJ(J)
SP=P(I,J)
DO 569 L=1,LS1-1
569 PLO(L)=SP*SIG(L)+PTOP
IF (PM(K).GE.SP+PTOP) GO TO 600
L=1
IF (PM(K).GE.PLO(1)) GO TO 580
570 LUP=L+1
IF (L.EQ.LM) GO TO 580
IF (PM(K).GE.PLO(LUP)) GO TO 575
L=LUP
GO TO 570
575 DPUP=PM(K)-PLO(LUP)
DPDN=PLO(L)-PM(K)
BYDP=1./(DPDN+DPUP)
TK=BYDP*(TX(I,J,L)*DPUP+TX(I,J,LUP)*DPDN)
QK=Q(I,J,L)*Q(I,J,LUP)/(BYDP*( * teeny)
ZK=BYDP*(PHI(I,J,L)*DPUP+PHI(I,J,LUP)*DPDN)
THK=BYDP*( GO TO 590
C**** SPECIAL CASES; L=1, L=LM
580 TK=TX(I,J,L)
QK=Q(I,J,L)
ZK=PHI(I,J,L)
THK=T(I,J,L)
C**** MERIDIONAL AVERAGING
590 WI=WI+W(I,J,K)
TKI=TKI+TK
QKI=QKI+QK
ZKI=ZKI+ZK
WTI=WTI+W(I,J,K)*TK
WQI=WQI+W(I,J,K)*QK
WZI=WZI+W(I,J,K)*ZK
THKI=THKI+THK
WTHI=WTHI+W(I,J,K)*THK
FIMI=FIMI+1.
600 CONTINUE
BYFIM=teeny
IF (FIMI.GT.teeny) BYFIM=1./FIMI
AJK(J,K-1,JK_ZMFVTDSE)=AJK(J,K-1,JK_ZMFVTDSE)+
& BYFIM*(SHA*TKI+ZKI)*WI
AJK(J,K-1,JK_TOTVTDSE)=AJK(J,K-1,JK_TOTVTDSE)+SHA*WTI+WZI
AJK(J,K-1,JK_ZMFVTLH)=AJK(J,K-1,JK_ZMFVTLH)+BYFIM*QKI*WI
AJK(J,K-1,JK_TOTVTLH)=AJK(J,K-1,JK_TOTVTLH)+WQI
AJK(J,K-1,JK_VTGEOEDDY)=AJK(J,K-1,JK_VTGEOEDDY)+WZI-BYFIM*WI*ZKI
C AJK(J,K-1,JK_BAREKEGEN)=AJK(J,K-1,JK_BAREKEGEN)+WTI-BYFIM*WI*TKI
WJK(J,K)=BYFIM*WI/DXYP(J)
WTJK(J,K)=BYFIM*(WTHI-BYFIM*WI*THKI)/DXYP(J)
AJK(J,K-1,JK_EDDVTPT)=AJK(J,K-1,JK_EDDVTPT)+WTJK(J,K)
610 CONTINUE
C****
C**** BAROCLINIC EDDY KINETIC ENERGY GENERATION
C****
DO 630 J=J_0,J_1
DO 630 K=1,KM
FIMI=0.
W2I=0.
PAI=0.
WPA2I=0.
DO 626 I=1,IMAXJ(J)
SP=P(I,J)
DO 611 L=1,LS1-1
611 PL(L)=SP*SIGE(L)+PTOP
PS=SP+PTOP
IF (PM(K+1).GE.PS) GO TO 626
L=1
PDN=PS
IF (PM(K).GE.PS) GO TO 614
PDN=PM(K)
612 IF (PM(K).GT.PL(L+1)) GO TO 614
L=L+1
GO TO 612
614 LUP=L
616 IF (PM(K+1).GE.PL(LUP+1)) GO TO 618
LUP=LUP+1
GO TO 616
618 CONTINUE
PTK=0.
C**** INTERPOLATE HERE
620 PUP=PL(L+1)
IF (LUP.EQ.L) PUP=PM(K+1)
DP=PDN-PUP
PTK=PTK+DP*TX(I,J,L)
IF (LUP.EQ.L) GO TO 622
L=L+1
PDN=PL(L)
GO TO 620
C**** ACCUMULATE HERE
622 FIMI=FIMI+1.
WUP=0.
IF (K.LT.KM) WUP=W(I,J,K+1)
W2I=W2I+W(I,J,K)+WUP
PY=PMO(K)
IF (PM(K).GE.PS) PY=.5*(PS+PM(K+1))
PAK=PTK/PY
PAI=PAI+PAK
WPA2I=WPA2I+( 626 CONTINUE
630 AJK(J,K,JK_BAREKEGEN)=AJK(J,K,JK_BAREKEGEN)-
& (WPA2I-W2I*PAI/(FIMI+teeny))
C****
C**** ACCUMULATE UV VERTICAL TRANSPORTS
C****
C**** DOUBLE POLAR WINDS
DO 640 K=1,KM
IF(GRID%HAVE_SOUTH_POLE) THEN
WSP=2.*W(1,1,K)/FIM
DO I=1,IM
W(I,1,K)=WSP
ENDDO
ENDIF
IF(GRID%HAVE_NORTH_POLE) THEN
WNP=2.*W(1,JM,K)/FIM
DO I=1,IM
W(I,JM,K)=WNP
ENDDO
ENDIF
640 CONTINUE
C P already halo'ed; no need CALL CHECKSUM(grid, P, __LINE__, __FILE__)
C P already halo'ed; no need CALL HALO_UPDATE(grid, P, FROM=SOUTH)
DO 710 J=J_0STG,J_1STG
UEARTH=RADIUS*OMEGA*COSV(J)
I=IM
DO 650 IP1=1,IM
PSEC(I)=.25*( 650 I=IP1
DO 710 K=2,KM
W4I=0.
UKI=0.
WU4I=0.
WKE4I=0.
FIMI=0.
I=IM
DO 700 IP1=1,IM
SP=PSEC(I)
DO 660 L=1,LS1-1
660 PLO(L)=SP*SIG(L)+PTOP
IF (PM(K).GE.SP+PTOP) GO TO 700
L=1
IF (PM(K).GE.PLO(1)) GO TO 680
670 LUP=L+1
IF (L.EQ.LM) GO TO 680
IF (PM(K).GE.PLO(LUP)) GO TO 675
L=LUP
GO TO 670
675 DPUP=PM(K)-PLO(LUP)
DPDN=PLO(L)-PM(K)
BYDP=1./(DPDN+DPUP)
UK=BYDP*( VK=BYDP*( GO TO 690
C**** SPECIAL CASES; L=1,L=LM
680 UK=U(I,J,L)
VK=V(I,J,L)
C**** MERIDIONAL AVERAGING
690 W4=W(I,J-1,K)+W(IP1,J-1,K)+W(I,J,K)+W(IP1,J,K)
W4I=W4I+W4
UKI=UKI+UK
WU4I=WU4I+W4*UK
WKE4I=WKE4I+W4*(UK*UK+VK*VK)
FIMI=FIMI+1.
700 I=IP1
BYFIM=1./(FIMI+teeny)
WUJK(J,K)=.25*(WU4I-W4I*UKI*BYFIM)*BYFIM/DXYV(J)
AJK(J,K-1,JK_TOTVTKE)=AJK(J,K-1,JK_TOTVTKE)+WKE4I
AJK(J,K-1,JK_VTAMEDDY)=AJK(J,K-1,JK_VTAMEDDY)+WU4I-BYFIM*W4I*UKI
710 AJK(J,K-1,JK_TOTVTAM)=AJK(J,K-1,JK_TOTVTAM)+WU4I !+W4I*UEARTH
C****
C**** POTENTIAL VORTICITY AND VERTICAL TRANSPORT OF POT. VORT.
C****
DO 760 J=J_0S,J_1S
JHEMI=1
IF (J.LT.1+JM/2) JHEMI=-1
DO 730 K=1,KM
PVI=0.
DO 720 I=1,IM
DUT(I,J,K)=JHEMI*STB(I,J,K)*(ZX(I,J,K)-FCOR(J))
720 PVI=PVI+DUT(I,J,K)
730 AJK(J,K,JK_POTVORT)=AJK(J,K,JK_POTVORT)+PVI
DO 760 K=2,KM
W2I=0.
PV2I=0.
WPV4I=0.
FIMI=0.
I=IM
DO 740 IP1=1,IM
PS=.5*( IF (PM(K).GE.PS) GO TO 740
W2=W(I,J,K)+W(IP1,J,K)
W2I=W2I+W2
PV2=DUT(I,J,K-1)+DUT(I,J,K)
PV2I=PV2I+PV2
WPV4I=WPV4I+W2*PV2
FIMI=FIMI+1.
740 I=IP1
AJK(J,K-1,JK_VTPV)=AJK(J,K-1,JK_VTPV)+WPV4I
760 AJK(J,K-1,JK_VTPVEDDY)=AJK(J,K-1,JK_VTPVEDDY)+
& WPV4I-W2I*PV2I/(FIMI+teeny)
C****
C**** SPECIAL MEAN/EDDY DIAGNOSTICS ARE CALCULATED
C****
DO 770 J=J_0STG,J_1STG
DO 765 K=2,KM
DPE=PMO(K)-PMO(K-1)
UP(J,K)=(UJK(J,K)-UJK(J,K-1))/DPE
765 PSIP(J,K)=(PSIJK(J,K)-PSIJK(J,K-1))/DPE
UP(J,1)=UP(J,2)
PSIP(J,1)=PSIP(J,2)
770 CONTINUE
DO 780 K=1,KM
KUP=K+1
IF (K.EQ.KM) KUP=KM
KDN=K-1
IF (K.EQ.1) KDN=1
CALL CHECKSUM(grid, TJK, 1614, "DIAG.f")
CALL HALO_UPDATE(grid, TJK, FROM=SOUTH)
DO 780 J=J_0STG,J_1STG
TY(J,K)=(TJK(J,K)-TJK(J-1,K))/DYV(J)
C**** E-P FLUX NORTHWARD COMPONENT
AJK(J,K,JK_EPFLXNCP)=AJK(J,K,JK_EPFLXNCP)+
& PSIJK(J,K)*(UJK(J,KUP)-UJK(J,KDN))/
* (PMO(KUP)-PMO(KDN))-UVJK(J,K)
780 CONTINUE
CALL CHECKSUM(grid, UJK, 1625, "DIAG.f")
CALL HALO_UPDATE(grid, UJK, FROM=NORTH)
CALL CHECKSUM(grid, PSIJK, 1627, "DIAG.f")
CALL HALO_UPDATE(grid, PSIJK, FROM=NORTH)
CALL CHECKSUM(grid, DXV, 1629, "DIAG.f")
CALL HALO_UPDATE(grid, DXV, FROM=NORTH)
CALL CHECKSUM(grid, VJK, 1631, "DIAG.f")
CALL HALO_UPDATE(grid, VJK, FROM=NORTH)
CALL CHECKSUM(grid, UP, 1633, "DIAG.f")
CALL HALO_UPDATE(grid, UP, FROM=NORTH)
CALL CHECKSUM(grid, TY, 1635, "DIAG.f")
CALL HALO_UPDATE(grid, TY, FROM=NORTH)
CALL CHECKSUM(grid, WJK, 1637, "DIAG.f")
CALL HALO_UPDATE(grid, WJK, FROM=NORTH)
CALL CHECKSUM(grid, PSIP, 1639, "DIAG.f")
CALL HALO_UPDATE(grid, PSIP, FROM=NORTH)
DO 800 J=J_0S,J_1S
DO 800 K=2,KM-1
UY=(UJK(J+1,K)*DXV(J+1)-UJK(J,K)*DXV(J)-FCOR(J))/DXYP(J)
PSIY=(PSIJK(J+1,K)*DXV(J+1)-PSIJK(J,K)*DXV(J))/DXYP(J)
C**** ZONAL MEAN MOMENTUM EQUATION (MEAN ADVECTION)
AJK(J,K,JK_DUDTMADV)=AJK(J,K,JK_DUDTMADV)-
& .5*UY*(VJK(J,K)+VJK(J+1,K))-
* .25*((UP(J+1,K+1)+UP(J,K+1))*WJK(J,K+1)+(UP(J+1,K)+UP(J,K))*
* WJK(J,K))
C**** ZONAL MEAN HEAT EQUATION (MEAN ADVECTION)
AJK(J,K,JK_DTDTMADV)=AJK(J,K,JK_DTDTMADV)-
& .5*(TY(J,K)*VJK(J,K)+TY(J+1,K)*VJK(J+1,K))
* -.5*STJK(J,K)*(WJK(J,K+1)+WJK(J,K))
C**** LAGRANGIAN MEAN MOMENTUM EQUATION (MEAN ADVECTION)
VSTAR=.5*(VJK(J,K)+VJK(J+1,K)-.5*(PSIP(J,K)+PSIP(J,K+1)
* +PSIP(J+1,K)+PSIP(J+1,K+1)))
WSTAR=.5*(WJK(J,K)+WJK(J,K+1))+PSIY
AJK(J,K,JK_DUDTTEM)=AJK(J,K,JK_DUDTTEM)-
& UY*VSTAR-.25*(UP(J,K)+UP(J+1,K)+
* UP(J,K+1)+UP(J+1,K+1))*WSTAR
AJK(J,K,JK_DTDTTEM)=AJK(J,K,JK_DTDTTEM)-
& .5*(TY(J+1,K)+TY(J,K))*VSTAR-
* STJK(J,K)*WSTAR
C**** VERTICAL E-P FLUX
AJK(J,K-1,JK_EPFLXVCP)=AJK(J,K-1,JK_EPFLXVCP)-
& WUJK(J,K)-.5*PSIJK(J,K)*UY
AJK(J,K,JK_EPFLXVCP)=AJK(J,K,JK_EPFLXVCP)-.5*PSIJK(J,K)*UY
800 CONTINUE
C****
C**** SPECTRAL ANALYSIS OF KINETIC ENERGIES AT CONSTANT PRESSURE
C****
IZERO=0
NM=1+IM/2
J45N=2+.75*(JM-1.)
c KS1=LS1
C**** TOTAL THE KINETIC ENERGIES
KE(:,:)=0.
C P already halo'ed; no need CALL CHECKSUM(grid, P, __LINE__, __FILE__)
C P already halo'ed; no need CALL HALO_UPDATE(grid, P, FROM=SOUTH)
DO 2140 J=J_0STG,J_1STG
I=IM
DO 2020 IP1=1,IM
PSEC(I)=.25*( 2020 I=IP1
DO 2140 K=1,KM
KSPHER=KLAYER(K)
IF (J.GT.JEQ) KSPHER=KSPHER+1
DO 2140 KX=IZERO,LM,LM
DO 2090 I=1,IM
DPUV=0.
SP=PSEC(I)
DO 2025 L=1,LS1-1
PLO(L)=SP*SIG(L)+PTOP ! PL or PLO ??
2025 PL(L)=SP*SIGE(L)+PTOP ! PLE or PL ??
PS=SP+PTOP
IF (PM(K+1).GE.PLO(1)) GO TO 2090 ! really ?? not PL?
L=1
PDN=PS
IF (PM(K).GE.PLO(1)) GO TO 2040 ! really ?? not PL?
PDN=PM(K)
2030 IF (PM(K).GT.PL(L+1)) GO TO 2040
L=L+1
GO TO 2030
2040 LUP=L
2050 IF (PM(K+1).GE.PL(LUP+1)) GO TO 2060
LUP=LUP+1
GO TO 2050
2060 CONTINUE
C**** ACCUMULATE HERE
SQRTDP=SQRT(PDN-PM(K+1))
2070 PUP=PL(L+1)
IF (LUP.EQ.L) PUP=PM(K+1)
DP=PDN-PUP
IF(KX.EQ.IZERO) DPUV=DPUV+DP*U(I,J,L)
IF(KX.EQ.LM) DPUV=DPUV+DP*V(I,J,L)
IF (LUP.EQ.L) GO TO 2080
L=L+1
PDN=PL(L)
GO TO 2070
2080 IF (SQRTDP.EQ.0.) SQRTDP=teeny
DPUV=DPUV/SQRTDP
2090 X1(I)=DPUV
CALL FFTE (X1,X1)
IF (J.EQ.JEQ) GO TO 2120
DO 2100 N=1,NM
2100 KE(N,KSPHER)=KE(N,KSPHER)+X1(N)*DXYV(J)
IF (J.NE.J45N) GO TO 2140
DO 2110 N=1,NM
2110 KE(N,KSPHER+2)=KE(N,KSPHER+2)+X1(N)*DXYV(J)
GO TO 2140
2120 DO 2130 N=1,NM
KE(N,KSPHER+2)=KE(N,KSPHER+2)+X1(N)*DXYV(J)
KE(N,KSPHER)=KE(N,KSPHER)+.5D0*X1(N)*DXYV(J)
2130 KE(N,KSPHER+1)=KE(N,KSPHER+1)+.5D0*X1(N)*DXYV(J)
2140 CONTINUE
DO 2150 KS=1,NSPHER
DO 2150 N=1,NM
2150 SPECA(N,18,KS)=SPECA(N,18,KS)+KE(N,KS)
C**** ACCUMULATE TIME USED IN DIAGA
CALL TIMEOUT(MBEGIN,MDIAG,MDYN)
RETURN
END SUBROUTINE DIAGB
SUBROUTINE DIAG7A 1,11
C****
C**** THIS ROUTINE ACCUMULATES A TIME SEQUENCE FOR SELECTED
C**** QUANTITIES AND FROM THAT PRINTS A TABLE OF WAVE FREQUENCIES.
C****
USE CONSTANT, only : grav,bygrav
USE MODEL_COM, only : im,imh,jm,lm,
& IDACC,JEQ,LS1,MDIAG,P,PTOP,PSFMPT,SIG,SIGE,U,V
USE DYNAMICS, only : PHI
USE DAGCOM, only : nwav_dag,wave,max12hr_sequ,j50n
USE DIAG_LOC, only : ldex
IMPLICIT NONE
REAL*8, DIMENSION(0:IMH) :: AN,BN
INTEGER, PARAMETER :: KM=6,KQMAX=12
INTEGER :: NMAX=nwav_dag
REAL*8, DIMENSION(IM,KM) :: HTRD
REAL*8, DIMENSION(KM), PARAMETER ::
& PMB=(/922.,700.,500.,300.,100.,10./),
& GHT=(/500.,2600.,5100.,8500.,15400.,30000./)
REAL*8 :: PIJ50N,PL,PLE,PLM1,SLOPE
INTEGER I,IDACC9,JLK,K,KQ,L,LX,MNOW,N
IDACC9=IDACC(9)+1
IDACC(9)=IDACC9
IF (IDACC9.GT.Max12HR_sequ) RETURN
DO KQ=1,3
CALL FFT ( DO N=1,NMAX
WAVE(1,IDACC9,N,2*KQ-1)=AN(N)
WAVE(2,IDACC9,N,2*KQ-1)=BN(N)
ENDDO
CALL FFT ( DO N=1,NMAX
WAVE(1,IDACC9,N,2*KQ)=AN(N)
WAVE(2,IDACC9,N,2*KQ)=BN(N)
ENDDO
ENDDO
DO 150 I=1,IM
PIJ50N=P(I,J50N)
K=1
L=1
PL=SIG(1)*P(I,J50N)+PTOP
130 L=L+1
IF(L.GE.LS1) PIJ50N=PSFMPT
PLM1=PL
PL=SIG(L)*PIJ50N+PTOP
IF (PMB(K).LT.PL.AND.L.LT.LM) GO TO 130
C**** ASSUME THAT PHI IS LINEAR IN LOG P
SLOPE=(PHI(I,J50N,L-1)-PHI(I,J50N,L))/LOG(PLM1/PL)
140 HTRD(I,K)=(PHI(I,J50N,L)+SLOPE*LOG(PMB(K)/PL))*BYGRAV-GHT(K)
IF (K.GE.KM) GO TO 150
K=K+1
IF (PMB(K).LT.PL.AND.L.LT.LM) GO TO 130
GO TO 140
150 CONTINUE
DO KQ=7,KQMAX
CALL FFT(HTRD(1,KQ-6),AN,BN)
DO N=1,NMAX
WAVE(1,IDACC9,N,KQ)=AN(N)
WAVE(2,IDACC9,N,KQ)=BN(N)
END DO
END DO
CALL TIMER (MNOW,MDIAG)
RETURN
END SUBROUTINE DIAG7A
SUBROUTINE DIAGCA (M) 13,22
!@sum DIAGCA Keeps track of the conservation properties of angular
!@+ momentum, kinetic energy, mass, total potential energy and water
!@auth Gary Russell/Gavin Schmidt
!@ver 1.0
USE MODEL_COM, only : mdiag,itime
USE DAGCOM, only : icon_AM,icon_KE,icon_MS,icon_TPE
* ,icon_WM,icon_LKM,icon_LKE,icon_EWM,icon_WTG,icon_HTG
* ,icon_OMSI,icon_OHSI,icon_OSSI,icon_LMSI,icon_LHSI,icon_MLI
* ,icon_HLI,title_con
USE SOIL_DRV, only: conserv_WTG,conserv_HTG
IMPLICIT NONE
!@var M index denoting from where DIAGCA is called
INTEGER, INTENT(IN) :: M
C****
C**** THE PARAMETER M INDICATES WHEN DIAGCA IS BEING CALLED
C**** M=1 INITIALIZE CURRENT QUANTITY
C**** 2 AFTER DYNAMICS
C**** 3 AFTER CONDENSATION
C**** 4 AFTER RADIATION
C**** 5 AFTER PRECIPITATION
C**** 6 AFTER LAND SURFACE (INCL. RIVER RUNOFF)
C**** 7 AFTER FULL SURFACE INTERACTION
C**** 8 AFTER FILTER
C**** 9 AFTER OCEAN DYNAMICS (from MAIN)
C**** 10 AFTER DAILY
C**** 11 AFTER OCEAN DYNAMICS (from ODYNAM)
C**** 12 AFTER OCEAN SUB-GRIDSCALE PHYS
C****
REAL*8, EXTERNAL :: conserv_AM,conserv_KE,conserv_MS,conserv_PE
* ,conserv_WM,conserv_EWM,conserv_LKM,conserv_LKE,conserv_OMSI
* ,conserv_OHSI,conserv_OSSI,conserv_LMSI,conserv_LHSI
* ,conserv_MLI,conserv_HLI
REAL*8 MNOW
INTEGER NT
C**** ATMOSPHERIC ANGULAR MOMENTUM
CALL conserv_DIAG(M,conserv_AM,icon_AM)
C**** ATMOSPHERIC KINETIC ENERGY
CALL conserv_DIAG(M,conserv_KE,icon_KE)
C**** ATMOSPHERIC MASS
CALL conserv_DIAG(M,conserv_MS,icon_MS)
C**** ATMOSPHERIC TOTAL POTENTIAL ENERGY
CALL conserv_DIAG(M,conserv_PE,icon_TPE)
C**** ATMOSPHERIC TOTAL WATER MASS
CALL conserv_DIAG(M,conserv_WM,icon_WM)
C**** ATMOSPHERIC TOTAL WATER ENERGY
CALL conserv_DIAG(M,conserv_EWM,icon_EWM)
C**** LAKE MASS AND ENERGY
CALL conserv_DIAG(M,conserv_LKM,icon_LKM)
CALL conserv_DIAG(M,conserv_LKE,icon_LKE)
C**** OCEAN ICE MASS, ENERGY, SALT
CALL conserv_DIAG(M,conserv_OMSI,icon_OMSI)
CALL conserv_DIAG(M,conserv_OHSI,icon_OHSI)
CALL conserv_DIAG(M,conserv_OSSI,icon_OSSI)
C**** LAKE ICE MASS, ENERGY
CALL conserv_DIAG(M,conserv_LMSI,icon_LMSI)
CALL conserv_DIAG(M,conserv_LHSI,icon_LHSI)
C**** GROUND WATER AND ENERGY
CALL conserv_DIAG(M,conserv_WTG,icon_WTG)
CALL conserv_DIAG(M,conserv_HTG,icon_HTG)
C**** LAND ICE MASS AND ENERGY
CALL conserv_DIAG(M,conserv_MLI,icon_MLI)
CALL conserv_DIAG(M,conserv_HLI,icon_HLI)
C**** OCEAN CALLS ARE DEALT WITH SEPARATELY
CALL DIAGCO (M)
C****
CALL TIMER (MNOW,MDIAG)
RETURN
END SUBROUTINE DIAGCA
module DIAG 4
contains
SUBROUTINE DIAGCD (M,UX,VX,DUT,DVT,DT1,PIT) 5,6
!@sum DIAGCD Keeps track of the conservation properties of angular
!@+ momentum and kinetic energy inside dynamics routines
!@auth Gary Russell
!@ver 1.0
USE CONSTANT, only : omega,mb2kg
USE MODEL_COM, only : im,jm,lm,fim,mdiag,mdyn
USE GEOM, only : cosv,radius,ravpn,ravps
USE DAGCOM, only : consrv
IMPLICIT NONE
C****
C**** THE PARAMETER M INDICATES WHEN DIAGCD IS BEING CALLED
C**** M=1 AFTER ADVECTION IN DYNAMICS
C**** 2 AFTER CORIOLIS FORCE IN DYNAMICS
C**** 3 AFTER PRESSURE GRADIENT FORCE IN DYNAMICS
C**** 4 AFTER STRATOS DRAG IN DYNAMICS
C**** 5 AFTER FLTRUV IN DYNAMICS
C**** 6 AFTER GRAVITY WAVE DRAG IN DYNAMICS
C****
!@var M index denoting from where DIAGCD is called
INTEGER, INTENT(IN) :: M
!@var DT1 current time step
REAL*8, INTENT(IN) :: DT1
!@var UX,VX current velocities
REAL*8, INTENT(IN), DIMENSION(IM,JM,LM) :: UX,VX
!@var DUT,DVT current momentum changes
REAL*8, INTENT(IN), DIMENSION(IM,JM,LM) :: DUT,DVT
!@var PIT current pressure tendency
REAL*8, INTENT(IN), OPTIONAL, DIMENSION(IM,JM) :: PIT
REAL*8, DIMENSION(JM) :: PI
INTEGER :: I,J,L,MBEGIN,N,IP1
LOGICAL dopit
REAL*8 :: DUTI,DUTIL,RKEI,RKEIL
INTEGER, DIMENSION(6) ::
* NAMOFM=(/2,3,4,5,6,7/), NKEOFM=(/14,15,16,17,18,19/)
CALL GETTIME(MBEGIN)
C****
C**** PRESSURE TENDENCY FOR CHANGE BY ADVECTION
C****
IF (M.eq.1) THEN
dopit=.true.
PI(1)=FIM*PIT(1,1)
PI(JM)=FIM*PIT(1,JM)
DO J=2,JM-1
PI(J)=SUM(PIT(:,J))
END DO
ELSE
PI=0.
dopit=.false.
END IF
C****
C**** CHANGE OF ANGULAR MOMENTUM AND KINETIC ENERGY BY VARIOUS
C**** PROCESSES IN DYNAMICS
C****
!$OMP PARALLEL DO PRIVATE (J,L,I,DUTIL,RKEIL,DUTI,RKEI,N)
DO J=2,JM
DUTIL=0.
RKEIL=0.
DO L=1,LM
DUTI=0.
RKEI=0.
DO I=1,IM
DUTI=DUTI+DUT(I,J,L)
RKEI=RKEI+(UX(I,J,L)*DUT(I,J,L)+VX(I,J,L)*DVT(I,J,L))
END DO
DUTIL=DUTIL+DUTI
RKEIL=RKEIL+RKEI
END DO
N=NAMOFM(M)
if (dopit) DUTIL=DUTIL+2.*DT1*RADIUS*OMEGA*COSV(J)*
* (PI(J-1)*RAVPN(J-1)+PI(J)*RAVPS(J))
CONSRV(J,N)=CONSRV(J,N)+DUTIL*COSV(J)*RADIUS*mb2kg
N=NKEOFM(M)
CONSRV(J,N)=CONSRV(J,N)+RKEIL*mb2kg
END DO
!$OMP END PARALLEL DO
C****
CALL TIMEOUT(MBEGIN,MDIAG,MDYN)
RETURN
END SUBROUTINE DIAGCD
end module DIAG
SUBROUTINE conserv_DIAG (M,CONSFN,ICON) 18,2
!@sum conserv_DIAG generic routine keeps track of conserved properties
!@auth Gary Russell/Gavin Schmidt
!@ver 1.0
USE MODEL_COM, only : jm
USE DAGCOM, only : consrv,nofm
IMPLICIT NONE
!@var M index denoting from where routine is called
INTEGER, INTENT(IN) :: M
!@var ICON index for the quantity concerned
INTEGER, INTENT(IN) :: ICON
!@var CONSFN external routine that calculates total conserved quantity
EXTERNAL CONSFN
!@var TOTAL amount of conserved quantity at this time
REAL*8, DIMENSION(JM) :: TOTAL
INTEGER :: I,J,NM,NI
C**** NOFM contains the indexes of the CONSRV array where each
C**** change is to be stored for each quantity. If NOFM(M,ICON)=0,
C**** no calculation is done.
C**** NOFM(1,ICON) is the index for the instantaneous value.
IF (NOFM(M,ICON).gt.0) THEN
C**** Calculate current value TOTAL
CALL CONSFN(TOTAL)
NM=NOFM(M,ICON)
NI=NOFM(1,ICON)
C**** Accumulate difference from last time in CONSRV(NM)
IF (M.GT.1) THEN
DO J=1,JM
CONSRV(J,NM)=CONSRV(J,NM)+(TOTAL(J)-CONSRV(J,NI))
END DO
END IF
C**** Save current value in CONSRV(NI)
DO J=1,JM
CONSRV(J,NI)=TOTAL(J)
END DO
END IF
RETURN
C****
END SUBROUTINE conserv_DIAG
SUBROUTINE conserv_AM(AM),3
!@sum conserv_AM calculates total atmospheric angular momentum
!@auth Gary Russell/Gavin Schmidt
!@ver 1.0
USE CONSTANT, only : omega,radius,mb2kg
USE MODEL_COM, only : im,jm,lm,fim,ls1,dsig,p,u,psfmpt,pstrat
USE GEOM, only : cosv,dxyn,dxys,dxyv
IMPLICIT NONE
REAL*8, DIMENSION(JM) :: AM,PI
INTEGER :: I,IP1,J,L
REAL*8 :: UMI,UMIL
C****
C**** ANGULAR MOMENTUM
C****
AM(1)=0.
PI(1)=FIM*P(1,1)
PI(JM)=FIM*P(1,JM)
DO J=2,JM-1
PI(J)=0.
DO I=1,IM
PI(J)=PI(J)+P(I,J)
END DO
END DO
DO J=2,JM
UMIL=0.
DO L=1,LM
UMI=0.
I=IM
DO IP1=1,IM
IF(L.LT.LS1) THEN
UMI=UMI+U(I,J,L)*((P(I,J-1)+P(IP1,J-1))*DXYN(J-1)
* +( ELSE
UMI=UMI+U(I,J,L)*(2.*PSFMPT*DXYV(J))
END IF
I=IP1
END DO
UMIL=UMIL+UMI*DSIG(L)
END DO
AM(J)=(RADIUS*OMEGA*COSV(J)*((PI(J-1)*DXYN(J-1)+PI(J)*DXYS(J))
* +FIM*PSTRAT*DXYV(J))+.5*UMIL)*COSV(J)*RADIUS*mb2kg
END DO
RETURN
C****
END SUBROUTINE conserv_AM
SUBROUTINE conserv_KE(RKE) 4,3
!@sum conserv_KE calculates total atmospheric kinetic energy
!@auth Gary Russell/Gavin Schmidt
!@ver 1.0
USE CONSTANT, only : mb2kg
USE MODEL_COM, only : im,jm,lm,fim,dsig,ls1,p,u,v,psfmpt
USE GEOM, only : dxyn,dxys,dxyv
IMPLICIT NONE
REAL*8, DIMENSION(JM) :: RKE
INTEGER :: I,IP1,J,L
REAL*8 :: RKEI,RKEIL
C****
C**** KINETIC ENERGY
C****
RKE(1)=0.
DO J=2,JM
RKEIL=0.
DO L=1,LM
RKEI=0.
I=IM
DO IP1=1,IM
IF(L.LT.LS1) THEN
RKEI=RKEI+( * *((P(I,J-1)+P(IP1,J-1))*DXYN(J-1)+( * ))*DXYS(J))
ELSE
RKEI=RKEI+( * (2.*PSFMPT*DXYV(J))
END IF
I=IP1
END DO
RKEIL=RKEIL+RKEI*DSIG(L)
END DO
RKE(J)=0.25*RKEIL*mb2kg
END DO
RETURN
C****
END SUBROUTINE conserv_KE
SUBROUTINE conserv_MS(RMASS),2
!@sum conserv_MA calculates total atmospheric mass
!@auth Gary Russell/Gavin Schmidt
!@ver 1.0
USE CONSTANT, only : mb2kg
USE MODEL_COM, only : im,jm,fim,p,pstrat
IMPLICIT NONE
REAL*8, DIMENSION(JM) :: RMASS
INTEGER :: I,J
C****
C**** MASS
C****
RMASS(1) =FIM*( RMASS(JM)=FIM*( DO J=2,JM-1
RMASS(J)=FIM*PSTRAT
DO I=1,IM
RMASS(J)=RMASS(J)+P(I,J)
END DO
RMASS(J)=RMASS(J)*mb2kg
END DO
RETURN
C****
END SUBROUTINE conserv_MS
SUBROUTINE conserv_PE(TPE) 4,4
!@sum conserv_TPE calculates total atmospheric potential energy
!@auth Gary Russell/Gavin Schmidt
!@ver 1.0
USE CONSTANT, only : sha,mb2kg
USE MODEL_COM, only : im,jm,lm,fim,t,p,ptop,zatmo
USE GEOM, only : imaxj
USE DYNAMICS, only : pk,pdsig
IMPLICIT NONE
REAL*8, DIMENSION(JM) :: TPE
INTEGER :: I,J,L
REAL*8 :: TPEI,TPEIL,SGEOI
C****
C**** TOTAL POTENTIAL ENERGY (J/m^2)
C****
DO J=1,JM
TPEIL=0.
DO L=1,LM
TPEI=0.
DO I=1,IMAXJ(J)
TPEI=TPEI+T(I,J,L)*PK(L,I,J)*PDSIG(L,I,J)
END DO
TPEIL=TPEIL+TPEI
END DO
SGEOI=0.
DO I=1,IMAXJ(J)
SGEOI=SGEOI+ZATMO(I,J)*( END DO
TPE(J)=(SGEOI+TPEIL*SHA)*mb2kg
END DO
TPE(1)=FIM*TPE(1)
TPE(JM)=FIM*TPE(JM)
RETURN
C****
END SUBROUTINE conserv_PE
SUBROUTINE conserv_WM(WATER),4
!@sum conserv_WM calculates total atmospheric water mass
!@auth Gary Russell/Gavin Schmidt
!@ver 1.0
USE CONSTANT, only : mb2kg
USE MODEL_COM, only : im,jm,lm,fim,wm,q
USE GEOM, only : imaxj
USE DYNAMICS, only : pdsig
IMPLICIT NONE
REAL*8, DIMENSION(JM) :: WATER
INTEGER :: I,J,L
C****
C**** TOTAL WATER MASS (kg/m^2)
C****
DO J=1,JM
WATER(J) = 0.
DO L=1,LM
DO I=1,IMAXJ(J)
WATER(J)=WATER(J)+( END DO
END DO
END DO
WATER(1) = FIM*WATER(1)
WATER(JM)= FIM*WATER(JM)
RETURN
C****
END SUBROUTINE conserv_WM
SUBROUTINE conserv_EWM(EWATER),5
!@sum conserv_EWM calculates total atmospheric water energy
!@auth Gary Russell/Gavin Schmidt
!@ver 1.0
USE CONSTANT, only : mb2kg,shv,grav,lhe
USE MODEL_COM, only : im,jm,lm,fim,wm,t,q,p
USE GEOM, only : imaxj
USE DYNAMICS, only : pdsig, pmid, pk
USE CLOUDS_COM, only : svlhx
IMPLICIT NONE
REAL*8, PARAMETER :: HSCALE = 7.8d0 ! km
REAL*8, DIMENSION(JM) :: EWATER
INTEGER :: I,J,L
REAL*8 W,EL
C****
C**** TOTAL WATER ENERGY (J/m^2)
C****
DO J=1,JM
EWATER(J) = 0.
DO L=1,LM
DO I=1,IMAXJ(J)
c this calculation currently only calculates latent heat
W =( EL=( * *mb2kg
EWATER(J)=EWATER(J)+EL !+W*(SHV*T(I,J,L)*PK(L,I,J)+GRAV
! * *HSCALE*LOG(P(I,J)/PMID(L,I,J)))
END DO
END DO
END DO
EWATER(1) = FIM*EWATER(1)
EWATER(JM)= FIM*EWATER(JM)
RETURN
C****
END SUBROUTINE conserv_EWM
SUBROUTINE DIAG5D (M5,NDT,DUT,DVT) 3,8
USE MODEL_COM, only : im,imh,jm,lm,fim,
& DSIG,JEQ,LS1,MDIAG,MDYN
USE DAGCOM, only : speca,nspher,klayer
USE DIAG_LOC