#include "rundeck_opts.h" !@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 : DIST_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(DIST_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) !hack hack hack! TX(:,:,:) = 0.d0 RETURN END SUBROUTINE ALLOC_DIAG_LOC SUBROUTINE DIAGA 2,49 !@sum DIAGA accumulate various diagnostics during dynamics !@auth Original Development Team !@ver 1.0 USE CONSTANT, only : grav,rgas,kapa,lhe,sha,bygrav,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,lm_req,req_fac_m,pmidl00 USE GEOM, only : areag,cosp,dlat,dxv,dxyn,dxyp,dxys,dxyv,dyp,fcor * ,imaxj,sinp,bydxyv,rapvn,rapvs USE RAD_COM, only : rqt USE DIAG_COM, only : aj=>aj_loc, aregj_loc, jreg, * apj=>apj_loc, ajl=>ajl_loc,asjl=>asjl_loc,ail,j50n,j70n,j5nuv * ,j5suv,j5s,j5n,aij=>aij_loc,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_qm,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,ij_p1000,ij_p925,ij_p700 * ,ij_p600,ij_p500 USE DYNAMICS, only : pk,phi,pmid,plij, pit,SD,pedn,am USE PBLCOM, only : tsavg USE DIAG_LOC, only : w,tx,lupa,ldna,jet,tjl0 USE DOMAIN_DECOMP, only : GET, CHECKSUM, HALO_UPDATE, & CHECKSUM_COLUMN, HALO_UPDATE_COLUMN, & GRID, SOUTH, NORTH, GLOBALSUM USE DOMAIN_DECOMP, only : AM_I_ROOT USE GETTIME_MOD IMPLICIT NONE REAL*8, DIMENSION(LM) :: GMEAN REAL*8, DIMENSION(grid%J_STRT_HALO:grid%J_STOP_HALO,LM) :: & GMEAN_part 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(IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: & PUV REAL*8, DIMENSION(LM_REQ) :: TRI REAL*8, DIMENSION(IM) :: THSEC,PSEC,SQRTP 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, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM) :: & THGM_part REAL*8, DIMENSION(LM):: THGM,PI0,AMI,DPI,PMI REAL*8, DIMENSION(LM+1):: PLEI REAL*8 :: & BBYGV,BYSDSG,DLNP,DLNP12,DLNP23,DBYSD, & DLNS,DP,DS,DT2,DTHDP,DU,DUDP,DUDX,DV,DXYPJ,ELX, * ESEPS,FPHI,GAMC,GAMM,GAMX,GMEANL,P4,P4I, & PDN,PE,PHIRI,PIBYIM,PIJ,PITIJ,PITMN, * PKE,PL,PRT,PU4I,PUV4I,PV4I,PVTHP, * ROSSX,SDMN,SDPU,SMALL,SP,SP2,SS,T4,THETA,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, SLP, PS, ZS REAL*8 :: gsum(IM, LM) REAL*8, dimension(:,:,:), allocatable :: TMP INTEGER :: J_0, J_1, J_0S, J_1S, J_0STG, J_1STG, J_0H LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE 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, & J_STRT_HALO=J_0H, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, & HAVE_NORTH_POLE=HAVE_NORTH_POLE) IDACC(4)=IDACC(4)+1 BYSDSG=1./(1.-SIGE(LM+1)) DLNP12=LOG(REQ_FAC_M(1)/REQ_FAC_M(2)) ! LOG(.75/.35) DLNP23=LOG(REQ_FAC_M(2)/REQ_FAC_M(3)) ! LOG(.35/.1) C**** C**** FILL IN HUMIDITY AND SIGMA DOT ARRAYS AT THE POLES C**** IF(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 ! HAVE_SOUTH_POLE IF(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 ! HAVE_NORTH_POLE C**** C**** CALCULATE PK AND TX, THE REAL TEMPERATURE C**** IF(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 ! HAVE_SOUTH_POLE IF(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 ! 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 HALO_UPDATE(grid, P, FROM=SOUTH) DO J=J_0STG,J_1STG I=IM DO IP1=1,IM PUV(I,J)=RAPVN(J-1)*(P(I,J-1)+P(IP1,J-1))+ * RAPVS( J)*(P(I, J)+P(IP1, 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 AREGJ_LOC(JR,J,J_TX1)=AREGJ_LOC(JR,J,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) PS=P(I,J)+PTOP ZS=BYGRAV*ZATMO(I,J) AIJ(I,J,IJ_SLP)=AIJ(I,J,IJ_SLP)+SLP(PS,TSAVG(I,J),ZS)-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)+(Q(I,J,L-1)-Q(I,J,L))*(PMB(K)-PL)/(PDN-PL) 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 C**** BEGIN AMIP IF((P(I,J)+PTOP).LT.1000.)AIJ(I,J,IJ_P1000)= * AIJ(I,J,IJ_P1000)+1. IF((P(I,J)+PTOP).LT.925.)AIJ(I,J,IJ_P925)=AIJ(I,J,IJ_P925)+1. IF((P(I,J)+PTOP).LT.700.)AIJ(I,J,IJ_P700)=AIJ(I,J,IJ_P700)+1. IF((P(I,J)+PTOP).LT.600.)AIJ(I,J,IJ_P600)=AIJ(I,J,IJ_P600)+1. IF((P(I,J)+PTOP).LT.500.)AIJ(I,J,IJ_P500)=AIJ(I,J,IJ_P500)+1. C**** END AMIP 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) AIJ(I,J,IJ_QM)=AIJ(I,J,IJ_QM)+Q(I,J,L)*AM(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)+(Q(I,J,L)+WM(I,J,L))*PIJ * *DSIG(L)*FTYPE(IT,I,J) END DO AREGJ_LOC(JR,J,J_QP)=AREGJ_LOC(JR,J,J_QP)+(Q(I,J,L)+WM(I,J,L * ))*PIJ*DSIG(L)*DXYPJ AREGJ_LOC(JR,J,J_TX)=AREGJ_LOC(JR,J,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 HALO_UPDATE(grid, TX, FROM=NORTH+SOUTH) 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=(T(I,J,LS1-1)-T(I,J,1))/(PHI(I,J,LS1-1)-PHI(I,J,1)+teeny) DO IT=1,NTYPE AJ(J,J_DTDGTR,IT)=AJ(J,J_DTDGTR,IT)+SS*FTYPE(IT,I,J) END DO AREGJ_LOC(JR,J,J_DTDGTR)=AREGJ_LOC(JR,J,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=(T(I,J,LSTR)-T(I,J,LS1-1))/((PHI(I,J,LSTR)-PHI(I,J,LS1-1)) * +teeny) DO IT=1,NTYPE AJ(J,J_DTSGST,IT)=AJ(J,J_DTSGST,IT)+SS*FTYPE(IT,I,J) END DO AREGJ_LOC(JR,J,J_DTSGST)=AREGJ_LOC(JR,J,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(pmidl00(lm)/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 HALO_UPDATE(grid, DUDVSQ, FROM=NORTH) DO J=J_0S,J_1S PIBYIM=PI(J)*BYIM call calc_vert_amp(PIBYIM,LS1-1,PI0,AMI,DPI,PLEI,PMI) DLNP=LOG(PMI(1)/PMI(LS1-1)) 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 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 HALO_UPDATE(grid, DUDVSQ, FROM=NORTH) DO J=J_0S,J_1S PIBYIM=PI(J)*BYIM DLNP=LOG((SIG(LS1-1)*PIBYIM+PTOP)/pmidl00(LSTR)) 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 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 HALO_UPDATE(grid, TIL, FROM=NORTH+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 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))* * (U(I,J,L)+U(I,J+1,L))*(Q(I,J,L)+Q(IP1,J,L))*DSIG(L)*.125 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, __LINE__, __FILE__) !Not necessary here, done above CALL HALO_UPDATE(grid, P, FROM=SOUTH) CALL HALO_UPDATE_COLUMN(grid, PLIJ, FROM=SOUTH) !Not necessary here, done above CALL CHECKSUM(grid, TX, __LINE__, __FILE__) !Not necessary here, done above CALL HALO_UPDATE(grid, TX, FROM=SOUTH) CALL HALO_UPDATE(grid, PHI, FROM=SOUTH) CALL HALO_UPDATE(grid, Q, 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*.25 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)*0.0625d0 SP2=PLIJ(L,IP1,J-1)+PLIJ(L,IP1,J) AIJ(IP1,J,IJ_PVQ)=AIJ(IP1,J,IJ_PVQ)+.125*SP2 * *(V(I,J,L)+V(IP1,J,L))*(Q(IP1,J-1,L)+Q(IP1,J,L))*DSIG(L) I=IP1 END DO AJL(J,L,JL_ZMFNTMOM)=AJL(J,L,JL_ZMFNTMOM)+.25*PU4I*PV4I/P4I AJL(J,L,JL_TOTNTMOM)=AJL(J,L,JL_TOTNTMOM)+.25*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(T(I,J,L+1),T(I,J,L)) 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_part(:,L)=0. DO 740 J=J_0,J_1 DO 720 I=1,IMAXJ(J) 720 SQRTP(I)=SQRT(P(I,J)) 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_part(J,L)=GMEAN_part(J,L)+ * (SIG(L)*P(I,J)+PTOP)*(T(I,J,LUP)-T(I,J,LDN))* * DXYP(J)/(P(I,J)*PK(L,I,J)) 740 CONTINUE C**** CALCULATE APE DO 760 L=1,LM IF(HAVE_SOUTH_POLE) THEN THJL(1,L)=THJL(1,L)*FIM THSQJL(1,L)=THSQJL(1,L)*FIM ENDIF IF(HAVE_NORTH_POLE) THEN THJL(JM,L)=THJL(JM,L)*FIM THSQJL(JM,L)=THSQJL(JM,L)*FIM ENDIF THGM_part(:,L)=0. DO J=J_0,J_1 THGM_part(J,L)=THJL(J,L)*DXYP(J) ENDDO 760 continue CALL GLOBALSUM(grid,THGM_part(:,1:LM),THGM(1:LM),ALL=.TRUE.) THGM=THGM/AREAG CALL GLOBALSUM(grid,GMEAN_part(:,1:LM),GMEAN(1:LM),ALL=.TRUE.) DO 761 L=1,LM LP1=LUPA(L) LM1=LDNA(L) GMEANL=GMEAN(L)/((SIG(LM1)-SIG(LP1))*AREAG) DO 761 J=J_0,J_1 761 AJL(J,L,JL_APE)=AJL(J,L,JL_APE)+ & (THSQJL(J,L)-2.*THJL(J,L)*THGM(L)+THGM(L)*THGM(L)* * 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 END DO CALL GLOBALSUM(grid, U(1:IM,:,1:LM), gsum, & jband=(/J5SUV,J5NUV/)) IF (AM_I_ROOT()) AIL(1:IM,1:LM,IL_UEQ) = & AIL(1:IM,1:LM,IL_UEQ)+gsum CALL GLOBALSUM(grid, V(1:IM,:,1:LM), gsum, & jband=(/J5SUV,J5NUV/)) IF (AM_I_ROOT()) AIL(1:IM,1:LM,IL_VEQ) = & AIL(1:IM,1:LM,IL_VEQ)+gsum CALL GLOBALSUM(grid, TX(1:IM,:,1:LM)-TF, gsum, & jband=(/J5S,J5N/)) IF (AM_I_ROOT()) AIL(1:IM,1:LM,IL_TEQ) = & AIL(1:IM,1:LM,IL_TEQ)+gsum allocate(TMP(1:IM, & grid%J_STRT_HALO:grid%J_STOP_HALO,1:LM)) DO L=1,LM DO J = MAX(J_0,J5S),MIN(J_1,J5N) DO I=1,IM TMP(I,J,L) = Q(I,J,L)/QSAT(TX(I,J,L),LHE,PMID(L,I,J)) END DO END DO END DO CALL GLOBALSUM(grid, TMP(1:IM,:,1:LM), gsum, & jband=(/J5S,J5N/)) deallocate(TMP) IF (AM_I_ROOT()) AIL(1:IM,1:LM,IL_QEQ) = & AIL(1:IM,1:LM,IL_QEQ)+gsum CALL GLOBALSUM(grid, TX(1:IM,:,1:LM)-TF, gsum, & jband=(/J50N,J50N/)) IF (AM_I_ROOT()) AIL(1:IM,1:LM,IL_T50N)= & AIL(1:IM,1:LM,IL_T50N)+gsum CALL GLOBALSUM(grid, U(1:IM,:,1:LM), gsum, & jband =(/J50N,J50N+1/)) IF (AM_I_ROOT()) AIL(1:IM,1:LM,IL_U50N)= & AIL(1:IM,1:LM,IL_U50N)+gsum CALL GLOBALSUM(grid, TX(1:IM,:,1:LM)-TF, gsum, & jband=(/J70N,J70N/)) IF (AM_I_ROOT()) AIL(1:IM,1:LM,IL_T70N)= & AIL(1:IM,1:LM,IL_T70N)+gsum CALL GLOBALSUM(grid, U(1:IM,:,1:LM), gsum, & jband =(/J70N,J70N+1/)) IF (AM_I_ROOT()) AIL(1:IM,1:LM,IL_U70N)= & AIL(1:IM,1:LM,IL_U70N)+gsum 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 END DO CALL GLOBALSUM(grid, W(1:IM,:,1:LM-1), gsum(1:IM,1:LM-1), & jband=(/J5S,J5N/)) IF (AM_I_ROOT()) AIL(1:IM,1:LM-1,IL_WEQ) = & AIL(1:IM,1:LM-1,IL_WEQ)+gsum(1:IM,1:LM-1) CALL GLOBALSUM(grid, W(1:IM,:,1:LM-1), gsum(1:IM,1:LM-1), & jband=(/J50N,J50N/)) IF (AM_I_ROOT()) AIL(1:IM,1:LM-1,IL_W50N) = & AIL(1:IM,1:LM-1,IL_W50N)+gsum(1:IM,1:LM-1) CALL GLOBALSUM(grid, W(1:IM,:,1:LM-1), gsum(1:IM,1:LM-1), & jband=(/J70N,J70N/)) IF (AM_I_ROOT()) AIL(1:IM,1:LM-1,IL_W70N) = & AIL(1:IM,1:LM-1,IL_W70N)+gsum(1:IM,1:LM-1) C**** C**** ELIASSEN PALM FLUX C**** C**** NORTHWARD TRANSPORT !Not necessary here, done above CALL CHECKSUM(grid, P, __LINE__, __FILE__) !Not necessary here, done above CALL HALO_UPDATE(grid, P, FROM=SOUTH) CALL HALO_UPDATE(grid, T, FROM=SOUTH) DO 868 J=J_0STG,J_1STG I=IM DO 862 IP1=1,IM PSEC(I)=(P(I,J )+P(IP1,J ))*RAPVS(J)+ * (P(I,J-1)+P(IP1,J-1))*RAPVN(J-1) 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, __LINE__, __FILE__) !Not necessary here, done above CALL HALO_UPDATE(grid, U, FROM=NORTH) 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)*(U(I,J+1,L)+U(I,J+1,L+1))-DXV(J)* * (U(I,J,L)+U(I,J,L+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*(T(I,J,L)+T(I,J,L+1)-THMN) 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)+.25* & ((.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 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 2,39 !@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,jdate, & mdyn,mdiag, ndaa,sig,sige,dsig,Jhour,u,v,t,p,q,wm USE GEOM, only : bydxyp,bydxyv,rapvs,rapvn, & COSV,DXV,DXYN,DXYP,DXYS,DXYV,DYP,DYV,FCOR,IMAXJ,RADIUS USE DIAG_COM, only : ajk=>ajk_loc,aijk=>aijk_loc,speca,nspher, ! adiurn,hdiurn & nwav_dag,ndiupt,hr_in_day,ijk_u,ijk_v,ijk_t,ijk_q,ijk_dp * ,ijk_dse,klayer,idd_w,ijdd,ijk_r,ijk_w,ijk_pf, * ijk_uv,ijk_vt,ijk_vq,ijk_vv,ijk_uu,ijk_tt, & 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,SD,pmid,pedn USE DIAG_LOC, only : w,tx,pm,pl,pmo,plo USE DOMAIN_DECOMP, only : GET, CHECKSUM, HALO_UPDATE, GRID USE DOMAIN_DECOMP, only : HALO_UPDATEj USE DOMAIN_DECOMP, only : SOUTH, NORTH, GLOBALSUM USE GETTIME_MOD IMPLICIT NONE REAL*8, DIMENSION(IMH+1,NSPHER) :: KE REAL*8, DIMENSION & (IMH+1,GRID%J_STRT_HALO:GRID%J_STOP_HALO,NSPHER) :: KE_part REAL*8, DIMENSION(IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM) :: & ZX,STB,UDX 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,P00,AML,PDSIGL,PMIDL REAL*8, DIMENSION(LM+1) :: PEDNL c REAL*8, DIMENSION(size(adiurn,1),size(adiurn,3), c & GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: ADIURN_part c REAL*8, DIMENSION(size(hdiurn,1),size(hdiurn,3), c & GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: HDIURN_part c REAL*8 :: ADIURNSUM,HDIURNSUM 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,KM REAL*8 :: & BYDP,BYFIM,DP,DPDN,DP4, & 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,SDK, & SMALL,SP,SQRTDP,THK,THKI,THPI,TK,TKI,TPI, & UDUTI, 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 & ,AMRHT,AMRHQ,AMUV,AMVQ,AMVT,AMUU,AMVV,AMTT REAL*8, PARAMETER :: BIG=1.E20 REAL*8 :: QSAT REAL*8 :: pm_ge_ps(im,grid%j_strt_halo:grid%j_stop_halo,lm) INTEGER :: J_0, J_1, J_0S, J_1S, J_0STG, J_1STG, J_0H LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE 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, & J_STRT_HALO=J_0H, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, & HAVE_NORTH_POLE=HAVE_NORTH_POLE) pm_ge_ps(:,:,:)=-1. C**** C**** INTERNAL QUANTITIES T,TH,Q,RH C**** KM=LM 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) call calc_vert_amp(SP,LM,P00,AML,PDSIGL,PEDNL,PMIDL) 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.PEDNL(L+1)) GO TO 120 L=L+1 GO TO 110 120 LUP=L 130 IF (PM(K+1).GE.PEDNL(LUP+1)) 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=PEDNL(L+1) 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,PMIDL(L)) 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=PEDNL(L) 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) SP=.5*(P(I,J)+P(IP1,J)) call calc_vert_amp(SP,LS1-1,P00,AML,PDSIGL,PEDNL,PMIDL) DO 175 L=1,LS1-1 PLO(L)=PMIDL(L) 175 PL(L)=PEDNL(L) DO 180 L=1,LM-1 DTH(L)=(T(I,J,L)+T(IP1,J,L)-T(I,J,L+1)-T(IP1,J,L+1))/ * (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*(P(I,J)+P(IP1,J))+PTOP 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**** ZX(:,:,:)=0. IF(HAVE_SOUTH_POLE) THEN UDX(:,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 HALO_UPDATE(grid, P, FROM=SOUTH) CALL HALO_UPDATE(grid, TX, FROM=SOUTH) CALL HALO_UPDATE(grid, PHI, FROM=SOUTH) CALL HALO_UPDATE(grid, Q, FROM=SOUTH) CALL HALO_UPDATE(grid, T, FROM=SOUTH) CALL HALO_UPDATEj(grid, STJK, FROM=SOUTH) c*** DO L=1,LM c*** CALL HALO_UPDATE(grid, STJK(:,L), FROM=SOUTH) c*** END DO DO 390 J=J_0STG,J_1STG I=IM DO 280 IP1=1,IM PSEC(I)=(P(I,J )+P(IP1,J ))*RAPVS(J)+ * (P(I,J-1)+P(IP1,J-1))*RAPVN(J-1) call calc_vert_amp(PSEC(I),LM,P00,AML,PDSIGL,PEDNL,PMIDL) DO K=1,KM UDX(I,J,K)=0. END DO DO L=1,LM DUT(I,J,L)=DUT(I,J,L)/(PDSIGL(L)*DXYV(J)) DVT(I,J,L)=DVT(I,J,L)/(PDSIGL(L)*DXYV(J)) END DO c DO 275 L=1,LS1-1 c DUT(I,J,L)=DUT(I,J,L)/(PSEC(I)*DXYV(J)*DSIG(L)) c 275 DVT(I,J,L)=DVT(I,J,L)/(PSEC(I)*DXYV(J)*DSIG(L)) c DO 276 L=LS1,LM c DUT(I,J,L)=DUT(I,J,L)/(PSFMPT*DXYV(J)*DSIG(L)) c 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) call calc_vert_amp(SP,LM,P00,AML,PDSIGL,PEDNL,PMIDL) PS=SP+PTOP DO 286 L=1,LS1-1 286 PL(L)=PEDNL(L) IF (PM(K+1).GE.PS) THEN pm_ge_ps(i,j,k) = 1. UDX(I,J,K)=BIG ELSE 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**** FOR AMIP AMVQ=0. AMVT=0. AMUU=0. AMVV=0. AMUV=0. AMTT=0. C**** END AMIP C**** INTERPOLATE HERE 330 PUP=PL(L+1) IF (LUP.EQ.L) PUP=PM(K+1) DP=PDN-PUP DP4=.25*DP PUK=PUK+DP*U(I,J,L) PVK=PVK+DP*V(I,J,L) PT4K=PT4K+(TX(I,J-1,L)+TX(IP1,J-1,L)+TX(I,J,L)+TX(IP1,J,L))*DP4 PZ4K=PZ4K+(PHI(I,J-1,L)+PHI(IP1,J-1,L)+PHI(I,J,L)+PHI(IP1,J,L)) & *DP4 PQ4K=PQ4K+(Q(I,J-1,L)+Q(IP1,J-1,L)+Q(I,J,L)+Q(IP1,J,L))*DP4 DUTK=DUTK+DP*DUT(I,J,L) DVTK=DVTK+DP*DVT(I,J,L) PS4K=PS4K+(T(I,J-1,L)+T(IP1,J-1,L)+T(I,J,L)+T(IP1,J,L))*DP4 C**** FOR AMIP 2 AMRHT=.25*(TX(I,J-1,L)+TX(IP1,J-1,L)+TX(I,J,L)+TX(IP1,J,L)) AMRHQ=.25*(Q(I,J-1,L)+Q(IP1,J-1,L)+Q(I,J,L)+Q(IP1,J,L)) AMVQ=AMVQ+DP*V(I,J,L)*AMRHQ AMVT=AMVT+DP*V(I,J,L)*AMRHT AMUU=AMUU+DP*U(I,J,L)*U(I,J,L) AMVV=AMVV+DP*V(I,J,L)*V(I,J,L) AMUV=AMUV+DP*U(I,J,L)*V(I,J,L) AMTT=AMTT+DP*AMRHT*AMRHT C**** END AMIP 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 AIJK(I,J,K,IJK_R) =AIJK(I,J,K,IJK_R) + * PQ4K/QSAT(PT4K/DPK,LHE,PMO(K)) AIJK(I,J,K,IJK_PF) =AIJK(I,J,K,IJK_PF)+1. C * * * FOR AMIP 2 * * * AIJK(I,J,K,IJK_UV)=AIJK(I,J,K,IJK_UV)+AMUV AIJK(I,J,K,IJK_VQ)=AIJK(I,J,K,IJK_VQ)+AMVQ AIJK(I,J,K,IJK_VT)=AIJK(I,J,K,IJK_VT)+AMVT AIJK(I,J,K,IJK_UU)=AIJK(I,J,K,IJK_UU)+AMUU AIJK(I,J,K,IJK_VV)=AIJK(I,J,K,IJK_VV)+AMVV AIJK(I,J,K,IJK_TT)=AIJK(I,J,K,IJK_TT)+AMTT C**** END AMIP C**** EDDY TRANSPORT OF THETA; VORTICITY 334 PS4I=PS4I+PS4K PSV4I=PSV4I+BYDP*PVK*PS4K UDX(I,J,K)=BYDP*PUK*DXV(J) !ESMF IF (UDX(I,J-1,K).LT.BIG) ZX(I,J-1,K)=UDX(I,J,K)-UDX(I,J-1,K) !ESMF IF (UDX(I,J-1,K).GE.BIG) ZX(I,J-1,K)=0. !ESMF IF (ZX(I,J-1,K).GE.BIG) ZX(I,J-1,K)=0. END IF !---> (PM(K+1).GE.PS) 340 I=IP1 !-->END I Loop (IP1) from IM to IM-1 (1-IM). 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)=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**** ZX for distributed parallelization c**** CALL HALO_UPDATE( grid, UDX, from=NORTH ) CALL HALO_UPDATE( grid, pm_ge_ps, from=NORTH) DO J=J_0,J_1S DO K=1,KM DO I=1,IM if (pm_ge_ps(i,j+1,k) < 0) then IF (UDX(I,J,K).LT.BIG ) ZX(I,J,K)=-UDX(I,J,K)+UDX(I,J+1,K) IF (UDX(I,J,K).GE.BIG) ZX(I,J,K)=0. IF (ZX(I,J,K).GE.BIG) ZX(I,J,K)=0 end if END DO END DO END DO C**** C**** alternate vertical mass flux diagnostic (from SD) C**** DO J=J_0,J_1 W(:,J,:)=0. END DO C**** interpolate SD to constant pressure DO J=J_0,J_1 I=IM DO IP1=1,IM DO K=1,KM-1 DPK=0. SDK=0. SP=P(I,J) DO L=1,LS1-1 PL(L)=PEDN(L,I,J) ! SP*SIGE(L)+PTOP END DO IF (PM(K+1).GE.SP+PTOP) GO TO 860 L=1 PDN=SP+PTOP IF (PM(K).GE.SP+PTOP) GO TO 820 PDN=PM(K) 810 IF (PM(K).GT.PL(L+1)) GO TO 820 L=L+1 GO TO 810 820 LUP=L 830 IF (PM(K+1).GE.PL(LUP+1)) GO TO 840 LUP=LUP+1 GO TO 830 840 CONTINUE C**** INTERPOLATE HERE 850 PUP=PL(L+1) IF (LUP.EQ.L) PUP=PM(K+1) DPK=DPK+(PDN-PUP) SDK=SDK+(PDN-PUP)*SD(I,J,L) IF (LUP.EQ.L) GO TO 860 L=L+1 PDN=PL(L) GO TO 850 860 CONTINUE C**** ACCUMULATE HERE (SHOULD I ACCUMULATE A WEIGHTING FUNCTION?) W(I,J,K)=0. IF (DPK.gt.0) THEN AIJK(I,J,K,IJK_W)=AIJK(I,J,K,IJK_W)+SDK*BYDXYP(J)/DPK W(I,J,K)=SDK/DPK END IF END DO I=IP1 END DO END DO 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 !!*** Warning: This diagnostic has 3 flaws (?) !!*** 1 - It assumes that DTsrc=1hr, (DTsrc=3600.) !!*** 2 - since DTdaa-Ndaa*DTsrc=2*DTdyn rather than 0, !!*** some hours are skipped once in a while !!*** 3 - Some of the first Ndaa hours are skipped at the !!*** beginning of a month and overcounted at the end; !!*** this happens to balance out, if and only if !!*** mod(days_in_month,ndaa)=0 (i.e. February if Ndaa=7) !!*** In addition, IHM occasionally is out-of-bounds. !! 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 I=1,IMAXJ(J) WI=WI+W(I,J,K) END DO 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)=PMID(L,I,J) ! 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*(Q(I,J,L)*DPDN+Q(I,J,LUP)*DPUP)+ * teeny) ZK=BYDP*(PHI(I,J,L)*DPUP+PHI(I,J,LUP)*DPDN) THK=BYDP*(T(I,J,L)*DPUP+T(I,J,LUP)*DPDN) 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)=PEDN(L,I,J) ! 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+(W(I,J,K)+WUP)*PAK 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(HAVE_SOUTH_POLE) THEN WSP=2.*W(1,1,K)/FIM DO I=1,IM W(I,1,K)=WSP ENDDO ENDIF IF(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) CALL HALO_UPDATE(grid, W, FROM=SOUTH) DO 710 J=J_0STG,J_1STG UEARTH=RADIUS*OMEGA*COSV(J) I=IM DO 650 IP1=1,IM PSEC(I)=.25*(P(I,J-1)+P(IP1,J-1)+P(I,J)+P(IP1,J)) 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*(U(I,J,L)*DPUP+U(I,J,LUP)*DPDN) VK=BYDP*(V(I,J,L)*DPUP+V(I,J,LUP)*DPDN) 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=.25*(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)=(WU4I-W4I*UKI*BYFIM)*BYFIM*BYDXYV(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*(P(I,J)+P(IP1,J))+PTOP IF (PM(K).GE.PS) GO TO 740 W2=.5*(W(I,J,K)+W(IP1,J,K)) W2I=W2I+W2 PV2=.5*(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 HALO_UPDATEj(grid, TJK, FROM=SOUTH) c*** DO L=1,LM c*** CALL HALO_UPDATE(grid, TJK(:,L), FROM=SOUTH) c*** END DO 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 HALO_UPDATEj(grid, PSIJK, FROM=NORTH) CALL HALO_UPDATEj(grid, UJK, FROM=NORTH) CALL HALO_UPDATEj(grid, VJK, FROM=NORTH) CALL HALO_UPDATEj(grid, WJK, FROM=NORTH) CALL HALO_UPDATEj(grid, UP, FROM=NORTH) CALL HALO_UPDATEj(grid, TY, FROM=NORTH) CALL HALO_UPDATEj(grid, PSIP, FROM=NORTH) c*** DO L=1,LM c*** CALL HALO_UPDATE(grid, PSIJK(:,L), FROM=NORTH) c*** CALL HALO_UPDATE(grid, UJK(:,L), FROM=NORTH) c*** CALL HALO_UPDATE(grid, VJK(:,L), FROM=NORTH) c*** If (L > 1) THEN c*** CALL HALO_UPDATE(grid, WJK(:,L), FROM=NORTH) c*** END IF c*** CALL HALO_UPDATE(grid, UP(:,L), FROM=NORTH) c*** CALL HALO_UPDATE(grid, TY(:,L), FROM=NORTH) c*** CALL HALO_UPDATE(grid, PSIP(:,L), FROM=NORTH) c*** END DO 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. KE_part(:,:,:)=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 J=J_0STG,J_1STG I=IM DO IP1=1,IM PSEC(I)=(P(I,J )+P(IP1,J ))*RAPVS(J)+ * (P(I,J-1)+P(IP1,J-1))*RAPVN(J-1) I=IP1 ENDDO DO K=1,KM KSPHER=KLAYER(K) IF (J.GT.JEQ) KSPHER=KSPHER+1 DO KX=IZERO,LM,LM DO I=1,IM DPUV=0. SP=PSEC(I) call calc_vert_amp(SP,LM,P00,AML,PDSIGL,PEDNL,PMIDL) DO 2025 L=1,LS1-1 PLO(L)=PMIDL(L) !SP*SIG(L)+PTOP ! PL or PLO ?? 2025 PL(L)=PEDNL(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 ENDDO CALL FFTE (X1,X1) IF (J.NE.JEQ) THEN DO N=1,NM KE_part(N,J,KSPHER)=KE_part(N,J,KSPHER)+X1(N)*DXYV(J) ENDDO IF (J.EQ.J45N) THEN DO N=1,NM KE_part(N,J,KSPHER+2)=KE_part(N,J,KSPHER+2)+ & X1(N)*DXYV(J) ENDDO ENDIF ELSE DO N=1,NM KE_part(N,J,KSPHER+2)=KE_part(N,J,KSPHER+2)+ & X1(N)*DXYV(J) KE_part(N,J,KSPHER )=KE_part(N,J,KSPHER)+ & .5D0*X1(N)*DXYV(J) KE_part(N,J,KSPHER+1)=KE_part(N,J,KSPHER+1)+ & .5D0*X1(N)*DXYV(J) ENDDO ENDIF ENDDO ENDDO ENDDO CALL GLOBALSUM(grid, KE_part, KE, ALL=.true.) 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,16 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,U,V USE DYNAMICS, only : PHI USE DIAG_COM, only : nwav_dag,wave,max12hr_sequ,j50n,kwp,re_and_im USE DIAG_LOC, only : ldex USE DOMAIN_DECOMP, only : GRID,GET,GLOBALSUM IMPLICIT NONE REAL*8, DIMENSION(0:IMH) :: AN,BN REAL*8, DIMENSION(RE_AND_IM,Max12HR_sequ,NWAV_DAG,KWP, & grid%J_STRT_HALO:grid%J_STOP_HALO) :: WAVE_part 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, DIMENSION(LM) :: P00,AML,PDSIGL,PMIDL REAL*8, DIMENSION(LM+1) :: PEDNL REAL*8 :: PIJ50N,PL,PLM1,SLOPE INTEGER I,IDACC9,K,KQ,L,MNOW,N INTEGER :: J_0, J_1 CALL GET(GRID,J_STRT=J_0,J_STOP=J_1) IDACC9=IDACC(9)+1 IDACC(9)=IDACC9 IF (IDACC9.GT.Max12HR_sequ) RETURN WAVE_part(:,:,:,1:6,:)=0. IF(J_0 <= JEQ .and. JEQ <= J_1) THEN DO KQ=1,3 CALL FFT (U(1,JEQ,LDEX(KQ)),AN,BN) DO N=1,NMAX WAVE_part(1,IDACC9,N,2*KQ-1,JEQ)=AN(N) WAVE_part(2,IDACC9,N,2*KQ-1,JEQ)=BN(N) ENDDO CALL FFT (V(1,JEQ,LDEX(KQ)),AN,BN) DO N=1,NMAX WAVE_part(1,IDACC9,N,2*KQ,JEQ)=AN(N) WAVE_part(2,IDACC9,N,2*KQ,JEQ)=BN(N) ENDDO ENDDO ENDIF CALL GLOBALSUM(grid, WAVE_part(1,IDACC9:IDACC9,1:NMAX,1:6,:), & WAVE(1,IDACC9:IDACC9,1:NMAX,1:6), ALL=.TRUE.) CALL GLOBALSUM(grid, WAVE_part(2,IDACC9:IDACC9,1:NMAX,1:6,:), & WAVE(2,IDACC9:IDACC9,1:NMAX,1:6), ALL=.TRUE.) WAVE_part(:,:,:,7:,:)=0. IF(J_0 <= J50N .and. J50N <= J_1) THEN DO 150 I=1,IM PIJ50N=P(I,J50N) call calc_vert_amp(PIJ50N,LM,P00,AML,PDSIGL,PEDNL,PMIDL) K=1 L=1 PL=PMIDL(1) ! SIG(1)*P(I,J50N)+PTOP 130 L=L+1 c IF(L.GE.LS1) PIJ50N=PSFMPT PLM1=PL PL=PMIDL(L) ! 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_part(1,IDACC9,N,KQ,J50N)=AN(N) WAVE_part(2,IDACC9,N,KQ,J50N)=BN(N) END DO END DO ENDIF CALL GLOBALSUM(grid, WAVE_part(1,IDACC9:IDACC9,1:NMAX,7:KQMAX,:), & WAVE(1,IDACC9:IDACC9,1:NMAX,7:KQMAX), ALL=.TRUE.) CALL GLOBALSUM(grid, WAVE_part(2,IDACC9:IDACC9,1:NMAX,7:KQMAX,:), & WAVE(2,IDACC9:IDACC9,1:NMAX,7:KQMAX), ALL=.TRUE.) CALL TIMER (MNOW,MDIAG) RETURN END SUBROUTINE DIAG7A SUBROUTINE DIAGCA (M) 17,24 !@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 #ifdef TRACERS_ON USE TRACER_COM, only: itime_tr0,ntm !xcon #endif USE DIAG_COM, 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**** 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 INTEGER 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) #ifdef TRACERS_ON C**** Tracer calls are dealt with separately do nt=1,ntm CALL DIAGTCA(M,NT) end do #endif C**** CALL TIMER (MNOW,MDIAG) RETURN END SUBROUTINE DIAGCA module DIAG 8 contains SUBROUTINE DIAGCD (grid,M,UX,VX,DUT,DVT,DT1,PIT) 9,11 !@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 DIAG_COM, only : consrv=>consrv_loc USE DOMAIN_DECOMP, only : GET, CHECKSUM, HALO_UPDATE, DIST_GRID USE DOMAIN_DECOMP, only : SOUTH USE GETTIME_MOD 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**** TYPE (DIST_GRID), INTENT(IN) :: grid !@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,GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM) :: & UX,VX !@var DUT,DVT current momentum changes REAL*8, INTENT(IN), & DIMENSION(IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM) :: & DUT,DVT !@var PIT current pressure tendency REAL*8, INTENT(IN), OPTIONAL, & DIMENSION(IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: PIT REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: 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/) INTEGER :: J_0, J_1, J_0S, J_1S, J_0STG, J_1STG, J_0H LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE 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, & J_STRT_HALO=J_0H, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, & HAVE_NORTH_POLE=HAVE_NORTH_POLE) C**** C**** PRESSURE TENDENCY FOR CHANGE BY ADVECTION C**** IF (M.eq.1) THEN dopit=.true. IF(HAVE_SOUTH_POLE) PI(1)=FIM*PIT(1,1) IF(HAVE_NORTH_POLE) PI(JM)=FIM*PIT(1,JM) DO J=J_0S,J_1S 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**** C**** CALL HALO_UPDATE(grid, PI, FROM=SOUTH) !$OMP PARALLEL DO PRIVATE (J,L,I,DUTIL,RKEIL,DUTI,RKEI,N) DO J=J_0STG,J_1STG 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) 23,4 !@sum conserv_DIAG generic routine keeps track of conserved properties !@auth Gary Russell/Gavin Schmidt !@ver 1.0 USE MODEL_COM, only : jm USE DOMAIN_DECOMP, only : GET, GRID, CHECKSUMj USE DIAG_COM, only : consrv=>consrv_loc,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(GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: TOTAL INTEGER :: I,J,NM,NI INTEGER :: J_0,J_1 CALL GET(grid, J_STRT=J_0, J_STOP=J_1) 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=J_0,J_1 CONSRV(J,NM)=CONSRV(J,NM)+(TOTAL(J)-CONSRV(J,NI)) END DO END IF C**** Save current value in CONSRV(NI) DO J=J_0,J_1 CONSRV(J,NI)=TOTAL(J) END DO END IF RETURN C**** END SUBROUTINE conserv_DIAG SUBROUTINE conserv_AM(AM),8 !@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 USE DOMAIN_DECOMP, only : GET, SOUTH, HALO_UPDATE, GRID USE DOMAIN_DECOMP, only : CHECKSUM IMPLICIT NONE REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: AM,PI INTEGER :: I,IP1,J,L REAL*8 :: UMI,UMIL INTEGER :: J_0S, J_1S, J_0STG, J_1STG LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE CALL GET(grid, J_STRT_SKP=J_0S, J_STOP_SKP=J_1S, & J_STRT_STGR=J_0STG, J_STOP_STGR=J_1STG, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, & HAVE_NORTH_POLE=HAVE_NORTH_POLE) C**** C**** ANGULAR MOMENTUM C**** IF(HAVE_SOUTH_POLE) PI(1)=FIM*P(1,1) IF(HAVE_SOUTH_POLE) AM(1)=0. IF(HAVE_NORTH_POLE) PI(JM)=FIM*P(1,JM) DO J=J_0S,J_1S PI(J)=0. DO I=1,IM PI(J)=PI(J)+P(I,J) END DO END DO CALL HALO_UPDATE(grid, PI, FROM=SOUTH) CALL HALO_UPDATE(grid, P, FROM=SOUTH) DO J=J_0STG,J_1STG 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) * +(P(I,J)+P(IP1,J))*DXYS(J)) 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) 1,7 !@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 USE DOMAIN_DECOMP, only : GET, CHECKSUM, HALO_UPDATE, GRID USE DOMAIN_DECOMP, only : SOUTH IMPLICIT NONE REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: RKE INTEGER :: I,IP1,J,L INTEGER :: J_0STG,J_1STG REAL*8 :: RKEI,RKEIL LOGICAL :: HAVE_SOUTH_POLE CALL GET(grid, J_STRT_STGR=J_0STG, J_STOP_STGR=J_1STG, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE) C**** C**** KINETIC ENERGY C**** CALL HALO_UPDATE(grid, P, FROM=SOUTH) IF (HAVE_SOUTH_POLE) RKE(1)=0. DO J=J_0STG,J_1STG RKEIL=0. DO L=1,LM RKEI=0. I=IM DO IP1=1,IM IF(L.LT.LS1) THEN RKEI=RKEI+(U(I,J,L)*U(I,J,L)+V(I,J,L)*V(I,J,L)) * *((P(I,J-1)+P(IP1,J-1))*DXYN(J-1)+(P(I,J)+P(IP1,J * ))*DXYS(J)) ELSE RKEI=RKEI+(U(I,J,L)*U(I,J,L)+V(I,J,L)*V(I,J,L))* * (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),4 !@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 USE DOMAIN_DECOMP, only : GET, GRID IMPLICIT NONE REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: RMASS INTEGER :: I,J INTEGER :: J_0S,J_1S LOGICAL :: HAVE_SOUTH_POLE,HAVE_NORTH_POLE CALL GET(grid, J_STRT_SKP=J_0S, J_STOP_SKP=J_1S, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, & HAVE_NORTH_POLE=HAVE_NORTH_POLE) C**** C**** MASS C**** IF(HAVE_SOUTH_POLE) RMASS(1) =FIM*(P(1,1) +PSTRAT)*mb2kg IF(HAVE_NORTH_POLE) RMASS(JM)=FIM*(P(1,JM)+PSTRAT)*mb2kg DO J=J_0S,J_1S 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) 1,6 !@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 USE DOMAIN_DECOMP, only : GET,GRID IMPLICIT NONE REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: TPE INTEGER :: I,J,L INTEGER :: J_0,J_1 LOGICAL :: HAVE_SOUTH_POLE,HAVE_NORTH_POLE REAL*8 :: TPEI,TPEIL,SGEOI CALL GET(grid, J_STRT=J_0, J_STOP=J_1, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, & HAVE_NORTH_POLE=HAVE_NORTH_POLE) C**** C**** TOTAL POTENTIAL ENERGY (J/m^2) C**** DO J=J_0,J_1 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)*(P(I,J)+PTOP) END DO TPE(J)=(SGEOI+TPEIL*SHA)*mb2kg END DO IF(HAVE_SOUTH_POLE) TPE(1)=FIM*TPE(1) IF(HAVE_NORTH_POLE) TPE(JM)=FIM*TPE(JM) RETURN C**** END SUBROUTINE conserv_PE SUBROUTINE conserv_WM(WATER),6 !@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 USE DOMAIN_DECOMP, only : GET, GRID IMPLICIT NONE REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: WATER INTEGER :: I,J,L INTEGER :: J_0,J_1 LOGICAL :: HAVE_NORTH_POLE, HAVE_SOUTH_POLE CALL GET(GRID, J_STRT=J_0, J_STOP=J_1, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, & HAVE_NORTH_POLE=HAVE_NORTH_POLE) C**** C**** TOTAL WATER MASS (kg/m^2) C**** DO J=J_0,J_1 WATER(J) = 0. DO L=1,LM DO I=1,IMAXJ(J) WATER(J)=WATER(J)+(Q(I,J,L)+WM(I,J,L))*PDSIG(L,I,J)*mb2kg END DO END DO END DO IF (HAVE_SOUTH_POLE) WATER(1) = FIM*WATER(1) IF (HAVE_NORTH_POLE) WATER(JM)= FIM*WATER(JM) RETURN C**** END SUBROUTINE conserv_WM SUBROUTINE conserv_EWM(EWATER),7 !@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 USE DOMAIN_DECOMP, only : GET, GRID IMPLICIT NONE REAL*8, PARAMETER :: HSCALE = 7.8d0 ! km REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: EWATER INTEGER :: I,J,L INTEGER :: J_0,J_1 LOGICAL :: HAVE_SOUTH_POLE,HAVE_NORTH_POLE REAL*8 W,EL CALL GET(GRID, J_STRT=J_0, J_STOP=J_1, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, & HAVE_NORTH_POLE=HAVE_NORTH_POLE) C**** C**** TOTAL WATER ENERGY (J/m^2) C**** DO J=J_0,J_1 EWATER(J) = 0. DO L=1,LM DO I=1,IMAXJ(J) c this calculation currently only calculates latent heat W =(Q(I,J,L)+WM(I,J,L))*PDSIG(L,I,J)*mb2kg EL=(Q(I,J,L)*LHE+WM(I,J,L)*(LHE-SVLHX(L,I,J)))*PDSIG(L,I,J) * *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 IF(HAVE_SOUTH_POLE) EWATER(1) = FIM*EWATER(1) IF(HAVE_NORTH_POLE) EWATER(JM)= FIM*EWATER(JM) RETURN C**** END SUBROUTINE conserv_EWM SUBROUTINE DIAG5D (M5,NDT,DUT,DVT) 5,13 USE MODEL_COM, only : im,imh,jm,lm,fim, & DSIG,JEQ,LS1,MDIAG,MDYN USE DIAG_COM, only : speca,nspher,klayer USE DIAG_LOC, only : FCUVA,FCUVB USE DOMAIN_DECOMP, only : GRID,GET,GLOBALSUM, WRITE_PARALLEL USE GETTIME_MOD IMPLICIT NONE REAL*8, DIMENSION(IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM) :: & DUT,DVT c REAL*8, DIMENSION(0:IMH,GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM,2) :: FCUVA,FCUVB c COMMON/WORK7/FCUVA,FCUVB INTEGER :: M5,NDT REAL*8, DIMENSION(IMH+1) :: X REAL*8, DIMENSION(0:IMH) :: FA,FB REAL*8, DIMENSION(IMH+1,NSPHER) :: KE REAL*8, DIMENSION & (IMH+1,GRID%J_STRT_HALO:GRID%J_STOP_HALO,NSPHER) :: KE_part INTEGER :: J,J45N,KUV,KSPHER,L,MBEGIN,MKE,N,NM INTEGER :: J_0STG,J_1STG CALL GET(GRID,J_STRT_STGR=J_0STG,J_STOP_STGR=J_1STG) NM=1+IM/2 J45N=2.+.75*(JM-1.) MKE=M5 GO TO (810,810,810,100,100, 100,810),M5 C**** 810 WRITE (6,910) M5 810 CALL WRITE_PARALLEL(M5, UNIT=6, format= & "('0INCORRECT VALUE OF M5 WHEN CALLING DIAG5D. M5=',I5)") C**** 910 FORMAT ('0INCORRECT VALUE OF M5 WHEN CALLING DIAG5D. M5=',I5) call stop_model('INCORRECT VALUE OF M5 WHEN CALLING DIAG5D',255) C**** C**** KINETIC ENERGY C**** C**** TRANSFER RATES FOR KINETIC ENERGY IN THE DYNAMICS 100 CALL GETTIME(MBEGIN) KE(:,:)=0. KE_part(:,:,:)=0. DO L=1,LM DO J=J_0STG,J_1STG KSPHER=KLAYER(L) IF (J > JEQ) KSPHER= KSPHER+1 DO KUV=1,2 ! loop over u,v IF(KUV.EQ.1) CALL FFT(DUT(1,J,L),FA,FB) IF(KUV.EQ.2) CALL FFT(DVT(1,J,L),FA,FB) DO N=1,NM X(N)=.5*FIM* & (FA(N-1)*FCUVA(N-1,J,L,KUV)+FB(N-1)*FCUVB(N-1,J,L,KUV)) ENDDO X(1)=X(1)+X(1) X(NM)=X(NM)+X(NM) IF (J.NE.JEQ) KE_part(:,J,KSPHER)=KE_part(:,J,KSPHER)+ & X(:)*DSIG(L) IF (J.EQ.J45N) THEN ! 45 N KE_part(:,J,KSPHER+2)=KE_part(:,J,KSPHER+2)+X(:)*DSIG(L) ELSE IF (J.EQ.JEQ) THEN ! EQUATOR DO N=1,NM KE_part(N,J,KSPHER+2)=KE_part(N,J,KSPHER+2)+ & X(N)*DSIG(L) KE_part(N,J,KSPHER )=KE_part(N,J,KSPHER )+ & .5D0*X(N)*DSIG(L) ! CONTRIB TO SH KE_part(N,J,KSPHER+1)=KE_part(N,J,KSPHER+1)+ & .5D0*X(N)*DSIG(L) ! CONTRIB TO NH ENDDO IF (KUV.EQ.2) KSPHER=KSPHER+1 ENDIF ENDDO ENDDO ENDDO CALL GLOBALSUM(grid, KE_part(1:NM,:,1:NSPHER), KE(1:NM,1:NSPHER), & ALL=.TRUE.) DO 180 KSPHER=1,NSPHER DO 180 N=1,NM 180 SPECA(N,MKE,KSPHER)=SPECA(N,MKE,KSPHER)+KE(N,KSPHER)/NDT C**** CALL TIMEOUT(MBEGIN,MDIAG,MDYN) RETURN END SUBROUTINE DIAG5D SUBROUTINE DIAG5A (M5,NDT) 10,20 C**** C**** THIS DIAGNOSTICS ROUTINE PRODUCES A SPECTRAL ANALYSIS OF KINETIC C**** AND AVAILABLE POTENTIAL ENERGIES AND THEIR TRANSFER RATES BY C**** VARIOUS ATMOSPHERIC PROCESSES. C**** C**** THE PARAMETER M INDICATES WHAT IS STORED IN SPECA(N,M,KSPHER), C**** IT ALSO INDICATES WHEN DIAG5A IS BEING CALLED. C**** M=1 MEAN STANDING KINETIC ENERGY BEFORE SOURCES C**** 2 MEAN KINETIC ENERGY BEFORE DYNAMICS C**** 3 MEAN POTENTIAL ENERGY C**** 4 CONVERSION OF K.E. BY ADVECTION AFTER ADVECTION C**** 5 CONVERSION OF K.E. BY CORIOLIS FORCE AFTER CORIOLIS TERM C**** 6 CONVERSION FROM P.E. INTO K.E. AFTER PRESS GRAD FORC C**** 7 CHANGE OF K.E. BY DYNAMICS AFTER DYNAMICS C**** 8 CHANGE OF P.E. BY DYNAMICS C**** 9 CHANGE OF K.E. BY CONDENSATION AFTER CONDENSATION C**** 10 CHANGE OF P.E. BY CONDENSATION C**** 11 CHANGE OF P.E. BY RADIATION AFTER RADIATION C**** 12 CHANGE OF K.E. BY SURFACE AFTER SURFACE C**** 13 CHANGE OF P.E. BY SURFACE C**** 14 CHANGE OF K.E. BY FILTER AFTER FILTER C**** 15 CHANGE OF P.E. BY FILTER C**** 16 CHANGE OF K.E. BY DAILY AFTER DAILY C**** 17 CHANGE OF P.E. BY DAILY C**** 18 UNUSED C**** 19 LAST KINETIC ENERGY C**** 20 LAST POTENTIAL ENERGY C**** USE CONSTANT, only : sha USE MODEL_COM, only : im,imh,jm,lm,fim, & DSIG,IDACC,JEQ,LS1,MDIAG, & P,PTOP,PSFMPT,SIG,T,U,V,ZATMO USE GEOM, only : AREAG,DXYN,DXYP,DXYS USE DIAG_COM, only : speca,atpe,nspher,kspeca,klayer USE DIAG_COM, only : SQRTM USE DYNAMICS, only : sqrtp,pk USE DOMAIN_DECOMP, only : GRID,GET,CHECKSUM,HALO_UPDATE, AM_I_ROOT USE DOMAIN_DECOMP, only : GLOBALSUM, SOUTH, WRITE_PARALLEL IMPLICIT NONE INTEGER :: M5,NDT REAL*8, DIMENSION(IM) :: X REAL*8, DIMENSION(IMH+1,NSPHER) :: KE,APE REAL*8, DIMENSION & (IMH+1,GRID%J_STRT_HALO:GRID%J_STOP_HALO,NSPHER) :: KE_part REAL*8, DIMENSION(IMH+1,4,LM) :: VAR REAL*8, DIMENSION(IMH+1,4,LM,GRID%J_STRT_HALO:GRID%J_STOP_HALO) & :: VAR_part REAL*8, DIMENSION(2) :: TPE REAL*8 :: TPE_sum REAL*8, DIMENSION(grid%J_STRT_HALO:grid%J_STOP_HALO) :: TPE_psum CMoved to DAGCOM so it could be declared allocatable REAL*8, SAVE, DIMENSION(IM,JM) :: SQRTM REAL*8, DIMENSION(LM) :: THJSP,THJNP,THGM INTEGER, PARAMETER :: IZERO=0 INTEGER, DIMENSION(KSPECA), PARAMETER :: & MTPEOF=(/0,0,1,0,0,0,0,2,0,3, 4,0,5,0,6,0,7,0,0,8/) INTEGER :: I,IJL2,IP1,J,J45N,JH,JHEMI,JP,K,KS,KSPHER,L,LDN, & LUP,MAPE,MKE,MNOW,MTPE,N,NM REAL*8 :: GMEAN(LM),GMTMP,SQRTPG,SUMI,SUMT,THGSUM,THJSUM REAL*8 :: THGSUM_part(grid%J_STRT_HALO:grid%J_STOP_HALO) REAL*8 :: GMSUM(grid%J_STRT_HALO:grid%J_STOP_HALO,LM) INTEGER :: J_0S,J_1S,J_0STG,J_1STG LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE CALL GET(GRID, J_STRT_SKP=J_0S , J_STOP_SKP=J_1S, & J_STRT_STGR=J_0STG, J_STOP_STGR=J_1STG, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, & HAVE_NORTH_POLE=HAVE_NORTH_POLE) SQRTPG = SQRT(PSFMPT) NM=1+IM/2 J45N=2.+.75*(JM-1.) IJL2=IM*JM*LM*2 MKE=M5 MAPE=M5 C**** C**** Note: KSPHER has been re-arranged from previous models to better C**** deal with optionally resolved stratosphere. The higher C**** values are only used if the model top is high enough. C**** C**** KSPHER=1 SOUTHERN TROPOSPHERE 2 NORTHERN TROPOSPHERE C**** 3 EQUATORIAL TROPOSPHERE 4 45 DEG NORTH TROPOSPHERE C**** C**** 5 SOUTHERN LOW STRATOSPHERE 6 NORTHERN LOW STRATOSPHERE C**** 7 EQUATORIAL LOW STRATOSPHERE 8 45 DEG NORTH LOW STRATOSPH C**** C**** 9 SOUTHERN MID STRATOSPHERE 10 NORTHERN MID STRATOSPHERE C**** 11 EQUATORIAL MID STRATOSPHERE 12 45 DEG NORTH MID STRATOSPH C**** C**** 13 SOUTHERN UPP STRATOSPHERE 14 NORTHERN UPP STRATOSPHERE C**** 15 EQUATORIAL UPP STRATOSPHERE 16 45 DEG NORTH UPP STRATOSPH C**** GO TO (200,200,810,810,810, 810,200,810,205,810, * 296,205,810,205,810, 205,810,810,810,810),M5 C*** 810 WRITE (6,910) M5 810 CALL WRITE_PARALLEL(M5, UNIT=6, format= & "('0INCORRECT VALUE OF M5 WHEN CALLING DIAG5A. M5=',I5)") C**** 910 FORMAT ('0INCORRECT VALUE OF M5 WHEN CALLING DIAG5A. M5=',I5) call stop_model('INCORRECT VALUE OF M5 WHEN CALLING DIAG5A.',255) C**** MASS FOR KINETIC ENERGY 200 CONTINUE I=IM CALL HALO_UPDATE(grid, P, FROM=SOUTH) DO J=J_0STG,J_1STG DO IP1=1,IM SQRTM(I,J)=SQRT(.5*((P(I,J)+P(IP1,J))*DXYS(J)+(P(I,J-1)+ * P(IP1,J-1))*DXYN(J-1))) I=IP1 END DO END DO C**** 205 CONTINUE MAPE=MKE+1 KE(:,:)=0. KE_part(:,:,:)=0. C**** CURRENT KINETIC ENERGY DO L=1,LM DO J=J_0STG,J_1STG IF (J <= JEQ) THEN KSPHER=KLAYER(L) ELSE KSPHER=KLAYER(L)+1 END IF DO K = IZERO,LM,LM IF(K.EQ.IZERO)X(1:IM)=U(1:IM,J,L)*SQRTM(1:IM,J) IF(K.EQ.LM) X(1:IM)=V(1:IM,J,L)*SQRTM(1:IM,J) CALL FFTE (X,X) IF (J.EQ.JEQ) THEN DO N=1,NM KE_part(N,J,KSPHER+2)=KE_part(N,J,KSPHER+2)+X(N)*DSIG(L) KE_part(N,J,KSPHER )=KE_part(N,J,KSPHER )+ & .5D0*X(N)*DSIG(L) KE_part(N,J,KSPHER+1)=KE_part(N,J,KSPHER+1)+ & .5D0*X(N)*DSIG(L) ENDDO cgsfc IF(K.EQ.LM)KSPHER=KSPHER+1 ELSE DO N=1,NM KE_part(N,J,KSPHER)=KE_part(N,J,KSPHER)+X(N)*DSIG(L) ENDDO IF (J.EQ.J45N) THEN DO N=1,NM KE_part(N,J,KSPHER+2)=KE_part(N,J,KSPHER+2)+ & X(N)*DSIG(L) ENDDO ENDIF ENDIF ENDDO ENDDO ENDDO CALL GLOBALSUM(grid, KE_part, KE, ALL=.true.) IF (NDT /= 0) THEN C**** TRANSFER RATES AS DIFFERENCES OF KINETIC ENERGY DO KS=1,NSPHER DO N=1,NM SPECA(N,MKE,KS)=SPECA(N,MKE,KS)+(KE(N,KS)-SPECA(N,19,KS))/NDT END DO END DO END IF DO KS=1,NSPHER DO N=1,NM SPECA(N,19,KS)=KE(N,KS) END DO END DO C**** C**** POTENTIAL ENERGY C**** 296 CONTINUE APE(:,:)=0. C**** CURRENT AVAILABLE POTENTIAL ENERGY DO L = 0, LM LDN = MAX(L - 1, 1) LUP = MIN(L + 1,LM) IF (L < LM) THEN IF(LUP.GE.LS1) THEN IF(HAVE_SOUTH_POLE) THJSP(LUP) = T(1,1,LUP)*SQRTPG IF(HAVE_NORTH_POLE) THJNP(LUP) = T(1,JM,LUP)*SQRTPG ELSE IF(HAVE_SOUTH_POLE) THJSP(LUP)=T(1,1,LUP)*SQRTP(1,1) IF(HAVE_NORTH_POLE) THJNP(LUP)=T(1,JM,LUP)*SQRTP(1,JM) ENDIF IF(HAVE_SOUTH_POLE) THGSUM_part(1) = FIM*THJSP(LUP)*DXYP(1) IF(HAVE_NORTH_POLE) THGSUM_part(JM) = FIM*THJNP(LUP)*DXYP(JM) DO J=J_0S,J_1S THJSUM=0. DO I=1,IM THJSUM=THJSUM+T(I,J,LUP)*SQRTP(I,J) ENDDO THGSUM_part(J) = THJSUM*DXYP(J) ENDDO CALL GLOBALSUM(grid, THGSUM_part, THGSUM, ALL=.TRUE.) THGM(LUP)=THGSUM/AREAG End IF IF (LUP < 2) CYCLE VAR(2:NM,1:2,L)=0. VAR_part(:,:,L,:)=0 IF(HAVE_SOUTH_POLE) THEN VAR_part(1,1,L,1)=.5*(THJSP(L)-THGM(L))**2*DXYP(1)*FIM GMSUM(1,L)=((THJSP(LUP)-THJSP(LDN))*DXYP(1)* * (SIG(L)*P(1,1)+PTOP)/ * (SQRTP(1,1)*P(1,1)*PK(L,1,1)))*FIM END IF IF(HAVE_NORTH_POLE) THEN VAR_part(1,2,L,JM)=.5*(THJNP(L)-THGM(L))**2*DXYP(JM)*FIM GMSUM(JM,L)=((THJNP(LUP)-THJNP(LDN))*DXYP(JM)* * (SIG(L)*P(1,JM)+PTOP)/(SQRTP(1,JM)*P(1,JM)*PK(L,1,JM)))*FIM END IF DO J=J_0S,J_1S IF (J < JEQ) THEN JHEMI = 1 ELSE JHEMI = 2 END IF GMTMP = 0 DO I=1,IM X(I)=T(I,J,L)*SQRTP(I,J)-THGM(L) GMTMP=GMTMP+(T(I,J,LUP)-T(I,J,LDN))*(SIG(L)*P(I,J)+PTOP)/ * (P(I,J)*PK(L,I,J)) END DO GMSUM(J,L) = GMTMP * DXYP(J) CALL FFTE (X,X) DO N=1,NM VAR_part(N,JHEMI,L,J)=VAR_part(N,JHEMI,L,J)+X(N)*DXYP(J) END DO IF (J == JEQ-1) THEN DO N=1,NM VAR_part(N,3,L,J)=X(N)*DXYP(J) END DO END IF IF (J == J45N-1) THEN DO N=1,NM VAR_part(N,4,L,J)=X(N)*DXYP(J) END DO END IF END DO END DO CALL GLOBALSUM(grid, GMSUM, GMEAN) CALL GLOBALSUM(grid, VAR_part, VAR) IF (AM_I_ROOT()) THEN DO L = 1, LM LDN = MAX(L - 1, 1) LUP = MIN(L + 1,LM) GMEAN(L)=DSIG(L)*AREAG*(SIG(LDN)-SIG(LUP))/GMEAN(L) KS=KLAYER(L) DO JHEMI=1,4 DO N=1,NM APE(N,KS)=APE(N,KS)+VAR(N,JHEMI,L)*GMEAN(L) END DO KS=KS+1 END DO END DO END IF C**** CURRENT TOTAL POTENTIAL ENERGY 450 CONTINUE IF (HAVE_SOUTH_POLE) THEN J=1 SUMT=0 DO L=1, LM SUMT=SUMT + T(1,J,L)*PK(L,1,J)*DSIG(L) END DO TPE_psum(J)=FIM*DXYP(J)*(ZATMO(1,J)*(P(1,J)+PTOP)+ * SUMT*SHA*P(1,J)) END IF IF (HAVE_NORTH_POLE) THEN J=JM SUMT=0 DO L=1, LM SUMT=SUMT + T(1,J,L)*PK(L,1,J)*DSIG(L) END DO TPE_psum(J)=FIM*DXYP(J)*(ZATMO(1,J)*(P(1,J)+PTOP)+ * SUMT*SHA*P(1,J)) END IF DO J=J_0S, J_1S SUMI=0 DO I=1,IM SUMT=0 DO L=1,LM SUMT=SUMT + T(I,J,L)*PK(L,I,J)*DSIG(L) END DO SUMI=SUMI+ZATMO(I,J)*(P(I,J)+PTOP)+SUMT*SHA*P(I,J) END DO TPE_psum(J) = SUMI*DXYP(J) END DO CALL GLOBALSUM(grid, TPE_psum, TPE_sum, TPE) IF (NDT /= 0) THEN MTPE=MTPEOF(MAPE) C**** TRANSFER RATES AS DIFFERENCES FOR POTENTIAL ENERGY DO KS=1,NSPHER DO N=1,NM SPECA(N,MAPE,KS)=SPECA(N,MAPE,KS)+(APE(N,KS)-SPECA(N,20,KS))/NDT END DO END DO ATPE(MTPE,1)=ATPE(MTPE,1)+(TPE(1)-ATPE(8,1))/NDT ATPE(MTPE,2)=ATPE(MTPE,2)+(TPE(2)-ATPE(8,2))/NDT END IF DO KS=1,NSPHER DO N=1,NM SPECA(N,20,KS)=APE(N,KS) END DO END DO ATPE(8,1)=TPE(1) ATPE(8,2)=TPE(2) IF (M5.EQ.2) THEN C**** ACCUMULATE MEAN KINETIC ENERGY AND MEAN POTENTIAL ENERGY DO KS=1,NSPHER DO N=1,NM SPECA(N,2,KS)=SPECA(N,2,KS)+KE(N,KS) SPECA(N,3,KS)=SPECA(N,3,KS)+APE(N,KS) END DO END DO ATPE(1,1)=ATPE(1,1)+TPE(1) ATPE(1,2)=ATPE(1,2)+TPE(2) END IF CALL TIMER (MNOW,MDIAG) RETURN END SUBROUTINE DIAG5A SUBROUTINE DIAG5F(UX,VX) 2,9 C**** FOURIER COEFFICIENTS FOR CURRENT WIND FIELD C**** USE MODEL_COM, only : im,imh,jm,lm, & IDACC,MDIAG,MDYN USE DIAG_LOC, only : FCUVA,FCUVB USE DOMAIN_DECOMP, only : GRID,GET USE GETTIME_MOD IMPLICIT NONE REAL*8, DIMENSION(IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM) :: & UX,VX c REAL*8, DIMENSION(0:IMH,JM,LM,2) :: FCUVA,FCUVB c COMMON/WORK7/FCUVA,FCUVB INTEGER :: J,L,MBEGIN INTEGER :: J_0STG, J_1STG CALL GET(GRID, J_STRT_STGR=J_0STG, J_STOP_STGR=J_1STG) CALL GETTIME(MBEGIN) IDACC(6)=IDACC(6)+1 DO L=1,LM DO J=J_0STG,J_1STG CALL FFT(UX(1,J,L),FCUVA(0,J,L,1),FCUVB(0,J,L,1)) CALL FFT(VX(1,J,L),FCUVA(0,J,L,2),FCUVB(0,J,L,2)) ENDDO ENDDO CALL TIMEOUT(MBEGIN,MDIAG,MDYN) RETURN END SUBROUTINE DIAG5F SUBROUTINE DIAG4A 1,2 C**** C**** THIS ROUTINE PRODUCES A TIME HISTORY OF ENERGIES C**** USE MODEL_COM, only : im,istrat,IDACC USE DIAG_COM, only : energy,speca,ned IMPLICIT NONE INTEGER :: I,IDACC5,N,NM IF (IDACC(4).LE.0.OR.IDACC(7).LE.0) RETURN NM=1+IM/2 C**** C**** LOAD ENERGIES INTO TIME HISTORY ARRAY C**** IDACC5=IDACC(5)+1 IF (IDACC5.GT.100) RETURN DO I=0,1+ISTRAT ! loop over number of 'spheres' ENERGY(1+NED*I,IDACC5)=SPECA(1,19,1+4*I) ! SH ENERGY(2+NED*I,IDACC5)=SPECA(1,19,2+4*I) ! NH ENERGY(5+NED*I,IDACC5)=SPECA(2,19,2+4*I) ! NH wave 1 ENERGY(6+NED*I,IDACC5)=SPECA(3,19,2+4*I) ! NH wave 2 ENERGY(7+NED*I,IDACC5)=SPECA(1,20,1+4*I) ENERGY(8+NED*I,IDACC5)=SPECA(1,20,2+4*I) DO N=2,NM ENERGY( 3+NED*I,IDACC5)=ENERGY( 3+10*I,IDACC5)+SPECA(N,19,1+4*I) ENERGY( 4+NED*I,IDACC5)=ENERGY( 4+10*I,IDACC5)+SPECA(N,19,2+4*I) ENERGY( 9+NED*I,IDACC5)=ENERGY( 9+10*I,IDACC5)+SPECA(N,20,1+4*I) ENERGY(10+NED*I,IDACC5)=ENERGY(10+10*I,IDACC5)+SPECA(N,20,2+4*I) END DO END DO IDACC(5)=IDACC5 RETURN C**** END SUBROUTINE DIAG4A module subdaily 1,6 !@sum SUBDAILY defines variables associated with the sub-daily diags !@auth Gavin Schmidt USE MODEL_COM, only : im,jm,lm,itime USE FILEMANAGER, only : openunit, closeunit, nameunit USE DIAG_COM, only : kgz_max,pmname,P_acc USE PARAM #ifdef TRACERS_ON USE TRACER_COM, only : ntm, trm, trname, mass2vol, n_Ox, n_SO4, * n_SO4_d1,n_SO4_d2,n_SO4_d3,n_clay,n_clayilli,n_sil1quhe, * n_water, n_HDO, n_Be7 #if (defined TRACERS_DUST) || (defined TRACERS_MINERALS) ||\ (defined TRACERS_QUARZHEM) * ,Ntm_dust #endif #ifdef TRACERS_DRYDEP & ,dodrydep #endif #ifdef TRACERS_WATER & ,dowetdep, trw0 #endif #ifdef TRACERS_COSMO USE COSMO_SOURCES, only : BE7D_acc,BE7W_acc #endif #endif IMPLICIT NONE SAVE !@var kddmax maximum number of sub-daily diags output files INTEGER, PARAMETER :: kddmax = 42 !@var kdd total number of sub-daily diags INTEGER :: kdd !@var kddunit total number of sub-daily files INTEGER :: kddunit !@var namedd array of names of sub-daily diags CHARACTER*10, DIMENSION(kddmax) :: namedd !@var iu_subdd array of unit numbers for sub-daily diags output INTEGER, DIMENSION(kddmax) :: iu_subdd !@var subddt = subdd + subdd1 = all variables for sub-daily diags CHARACTER*128 :: subddt = " " !@dbparam subdd string contains variables to save for sub-daily diags !@dbparam subdd1 additional string of variables for sub-daily diags C**** Note: for longer string increase MAX_CHAR_LENGTH in PARAM CHARACTER*64 :: subdd = "SLP" CHARACTER*64 :: subdd1 = " " !@dbparam Nsubdd: DT_save_SUBDD = Nsubdd*DTsrc sub-daily diag freq. INTEGER :: Nsubdd = 0 !@dbparam LmaxSUBDD: the max L when writing "ALL" levels INTEGER :: LmaxSUBDD = LM contains subroutine init_subdd(aDATE) 1,7 !@sum init_subdd initialise sub daily diags and position files !@auth Gavin Schmidt implicit none character*14, intent(in) :: adate character*12 name integer :: i,j,k,kunit,kk call sync_param( "subdd" ,subdd) call sync_param( "subdd1" ,subdd1) call sync_param( "Nsubdd",Nsubdd) call sync_param( "LmaxSUBDD",LmaxSUBDD) if (nsubdd.ne.0) then C**** combine strings subdd1 and subdd2: subddt=trim(subdd)//' '//subdd1 C**** calculate how many names k=0 i=1 10 j=index(subddt(i:len(subddt))," ") if (j.gt.1) then k=k+1 i=i+j else i=i+1 end if if (i.lt.len(subddt)) goto 10 kdd=k if (kdd.gt.kddmax) call stop_model * ("Increase kddmax: No. of sub-daily diags too big",255) C**** make array of names read(subddt,*) namedd(1:kdd) C**** open units and position call open_subdd(aDATE) C**** position correctly do kunit=1,kddunit call io_POS(iu_SUBDD(kunit),Itime-1,im*jm,Nsubdd) end do end if C**** initialise special subdd accumulation P_acc=0. #ifdef TRACERS_COSMO BE7W_acc=0. BE7D_acc=0. #endif return end subroutine init_subdd subroutine open_subdd(aDATE) 2,3 !@sum open_subdd opens sub daily diag files !@auth Gavin Schmidt implicit none character*14, intent(in) :: adate character*12 name integer :: k,kunit,kk kunit=0 do k=1,kdd C**** Some names have more than one unit associated (i.e. "ZALL") if (namedd(k)(len_trim(namedd(k))-2:len_trim(namedd(k))).eq. * "ALL") then select case (namedd(k)(1:1)) case ("U", "V", "W", "C", "D", "O", "B") ! velocities/tracers on model layers kunit=kunit+1 write(name,'(A1,A3,A7)') namedd(k)(1:1),'ALL',aDATE(1:7) call openunit(name,iu_SUBDD(kunit),.true.,.false.) case ("Z", "T", "R", "Q") ! heights, temps, rel/spec hum PMB levels do kk=1,kgz_max kunit=kunit+1 call openunit(namedd(k)(1:1)//trim(PMNAME(kk))// * aDATE(1:7),iu_SUBDD(kunit),.true.,.false.) end do end select else ! single file per name kunit=kunit+1 call openunit(trim(namedd(k))//aDATE(1:7),iu_SUBDD(kunit), * .true.,.false.) endif end do kddunit=kunit C**** return end subroutine open_subdd subroutine reset_subdd(aDATE) 1,2 !@sum reset_subdd resets sub daily diag files !@auth Gavin Schmidt implicit none character*14, intent(in) :: adate character*12 name if (nsubdd.ne.0) then C**** close and re-open units call closeunit ( iu_SUBDD(1:kddunit) ) call open_subdd( aDATE ) end if C**** return end subroutine reset_subdd subroutine get_subdd 1,94 !@sum get_SUBDD saves instantaneous variables at sub-daily frequency !@+ every ABS(NSUBDD) !@+ Note that TMIN, TMAX can only be output once a day !@+ Current options: SLP, PS, SAT, PREC, QS, LCLD, MCLD, HCLD, PTRO !@+ QLAT, QSEN, SWD, SWU, LWD, LWU, LWT, STX, STY, !@+ ICEF, SNOWD, TCLD, SST, SIT, US, VS, TMIN, TMAX !@+ Z*, R*, T*, Q* (on any fixed pressure level) !@+ U*, V*, W*, C* (on any model level) !@+ Ox* (on any model level with chemistry) !@+ D* (HDO on any model level) !@+ B* (BE7 on any model level) !@+ SO4 !@+ 7BEW, 7BED, BE7ATM, EVAP !@+ CTEM,CD3D,CI3D,CL3D,CDN3D,CRE3D,CLWP ! aerosol !@+ TAUSS,TAUMC,CLDSS,CLDMC !@+ SO4_d1,SO4_d2,SO4_d3, ! het. chem !@+ Clay, Silt1, Silt2, Silt3 ! dust !@+ DUEMIS,DUDEPTURB,DUDEPGRAV,DUDEPWET,DUTRS,DULOAD !@+ DUEMIS2 !@+ !@+ More options can be added as extra cases in this routine !@auth Gavin Schmidt/Reto Ruedy USE CONSTANT, only : grav,rgas,bygrav,bbyg,gbyrb,sday,tf,mair,sha * ,lhe,rhow,undef,stbo USE MODEL_COM, only : lm,p,ptop,zatmo,dtsrc,u,v,focean * ,flice,nday USE GEOM, only : imaxj,dxyp,bydxyp USE PBLCOM, only : tsavg,qsavg,usavg,vsavg USE CLOUDS_COM, only : llow,lmid,lhi,cldss,cldmc,taumc,tauss,fss #ifdef CLD_AER_CDNC * ,ctem,cd3d,ci3d,cl3d,cdn3d,cre3d,clwp #endif USE DYNAMICS, only : ptropo,am,wsave USE FLUXES, only : prec,dmua,dmva,tflux1,qflux1,uflux1,vflux1 * ,gtemp #if (defined TRACERS_DUST) || (defined TRACERS_MINERALS) ||\ (defined TRACERS_QUARZHEM) & ,dust_flux_glob,trs_glob #ifdef TRACERS_DRYDEP & ,depo_turb_glob,depo_grav_glob #endif #ifdef TRACERS_WATER & ,trprec #else & ,trprec_dust #endif #endif #ifdef TRACERS_DUST & ,dust_flux2_glob #endif USE SEAICE_COM, only : rsi,snowi USE LANDICE_COM, only : snowli USE LAKES_COM, only : flake USE GHY_COM, only : snowe,fearth USE RAD_COM, only : trhr,srdn,salb,cfrac,cosz1 USE DIAG_COM, only : z_inst,rh_inst,t_inst,kgz_max,pmname,tdiurn * ,p_acc,pmb USE DOMAIN_DECOMP, only : GRID,GET,WRITEI_PARALLEL IMPLICIT NONE REAL*4, DIMENSION(IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: DATA INTEGER :: I,J,K,L,kp,kunit,n,n1,n_fidx CHARACTER namel*3 REAL*8 POICE,PEARTH,PLANDI,POCEAN,QSAT,PS,SLP, ZS INTEGER :: J_0,J_1 LOGICAL :: HAVE_SOUTH_POLE,HAVE_NORTH_POLE CALL GET(GRID,J_STRT=J_0,J_STOP=J_1, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, & HAVE_NORTH_POLE=HAVE_NORTH_POLE) #ifdef TRACERS_DUST n_fidx=n_clay #else #ifdef TRACERS_MINERALS n_fidx=n_clayilli #else #ifdef TRACERS_QUARZHEM n_fidx=n_sil1quhe #endif #endif #endif kunit=0 C**** depending on namedd string choose what variables to output do k=1,kdd C**** simple diags select case (namedd(k)) case ("SLP") ! sea level pressure (mb) do j=J_0,J_1 do i=1,imaxj(j) ps=(p(i,j)+ptop) zs=bygrav*zatmo(i,j) data(i,j)=slp(ps,tsavg(i,j),zs) end do end do case ("PS") ! surface pressure (mb) do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=p(i,j)+ptop end do end do case ("SAT") ! surf. air temp (C) data=tsavg-tf case ("US") ! surf. u wind (m/s) data=usavg case ("VS") ! surf. v wind (m/s) data=vsavg case ("SST") ! sea surface temp (C) do j=J_0,J_1 do i=1,imaxj(j) if (FOCEAN(I,J)+FLAKE(I,J).gt.0) then data(i,j)=GTEMP(1,1,i,j) else data(i,j)=undef end if end do end do case ("SIT") ! surface ice temp (C) do j=J_0,J_1 do i=1,imaxj(j) if (RSI(I,J)*(FOCEAN(I,J)+FLAKE(I,J)).gt.0) then data(i,j)=GTEMP(1,2,i,j) else data(i,j)=undef end if end do end do case ("QS") ! surf humidity (kg/kg) data=qsavg case ("PREC") ! precip (mm/day) c data=sday*prec/dtsrc data=sday*P_acc/(Nsubdd*dtsrc) ! accum over Nsubdd steps P_acc=0. case ("SNOWD") ! snow depth (w.e. mm) do j=J_0,J_1 do i=1,imaxj(j) POICE=RSI(I,J)*(FOCEAN(I,J)+FLAKE(I,J)) PEARTH=FEARTH(I,J) PLANDI=FLICE(I,J) data(i,j)=1d3*(SNOWI(I,J)*POICE+SNOWLI(I,J)*PLANDI+SNOWE(I * ,J)*PEARTH)/RHOW end do end do case ("QLAT") ! latent heat (W/m^2) data=qflux1*lhe case ("QSEN") ! sensible heat flux (W/m^2) data=tflux1*sha case ("SWD") ! solar downward flux at surface (W/m^2) data=srdn*cosz1 ! multiply by instant cos zenith angle case ("SWU") ! solar upward flux at surface (W/m^2) ! estimating this from the downward x albedo, since that's already saved data=srdn*(1.-salb)*cosz1 case ("LWD") ! LW downward flux at surface (W/m^2) data=TRHR(0,:,:) case ("LWU") ! LW upward flux at surface (W/m^2) do j=J_0,J_1 do i=1,imaxj(j) POCEAN=(1.-RSI(I,J))*(FOCEAN(I,J)+FLAKE(I,J)) POICE=RSI(I,J)*(FOCEAN(I,J)+FLAKE(I,J)) PEARTH=FEARTH(I,J) PLANDI=FLICE(I,J) data(i,j)=STBO*(POCEAN*(GTEMP(1,1,I,J)+TF)**4+ * POICE *(GTEMP(1,2,I,J)+TF)**4+ * PLANDI*(GTEMP(1,3,I,J)+TF)**4+ * PEARTH*(GTEMP(1,4,I,J)+TF)**4) end do end do case ("LWT") ! LW upward flux at TOA (P1) (W/m^2) do j=J_0,J_1 ! sum up all cooling rates + net surface emission do i=1,imaxj(j) POCEAN=(1.-RSI(I,J))*(FOCEAN(I,J)+FLAKE(I,J)) POICE=RSI(I,J)*(FOCEAN(I,J)+FLAKE(I,J)) PEARTH=FEARTH(I,J) PLANDI=FLICE(I,J) data(i,j)=-SUM(TRHR(0:LM,I,J))+ * STBO*(POCEAN*(GTEMP(1,1,I,J)+TF)**4+ * POICE *(GTEMP(1,2,I,J)+TF)**4+ * PLANDI*(GTEMP(1,3,I,J)+TF)**4+ * PEARTH*(GTEMP(1,4,I,J)+TF)**4) end do end do case ("ICEF") ! ice fraction over open water (%) data=RSI*100. case ("STX") ! E-W surface stress (N/m^2) data=uflux1 case ("STY") ! N-S surface stress (N/m^2) data=vflux1 case ("LCLD") ! low level cloud cover (%) data=0. ! Warning: these can be greater >100! do j=J_0,J_1 do i=1,imaxj(j) do l=1,llow data(i,j)=data(i,j)+(cldss(l,i,j)+cldmc(l,i,j)) end do data(i,j)=data(i,j)*100./real(llow,kind=8) end do end do case ("MCLD") ! mid level cloud cover (%) data=0. ! Warning: these can be greater >100! do j=J_0,J_1 do i=1,imaxj(j) do l=llow+1,lmid data(i,j)=data(i,j)+(cldss(l,i,j)+cldmc(l,i,j)) end do data(i,j)=data(i,j)*100./real(lmid-llow,kind=8) end do end do case ("HCLD") ! high level cloud cover (%) data=0. ! Warning: these can be greater >100! do j=J_0,J_1 do i=1,imaxj(j) do l=lmid+1,lhi data(i,j)=data(i,j)+(cldss(l,i,j)+cldmc(l,i,j)) end do data(i,j)=data(i,j)*100./real(lhi-lmid,kind=8) end do end do case ("TCLD") ! total cloud cover (%) (As seen by rad) data=cfrac*100. case ("PTRO") ! tropopause pressure (mb) data = ptropo case ("TMIN") ! min daily temp (C) if (mod(itime+1,Nday).ne.0) then ! only at end of day kunit=kunit+1 cycle end if data=tdiurn(:,:,9) case ("TMAX") ! max daily temp (C) if (mod(itime+1,Nday).ne.0) then ! only at end of day kunit=kunit+1 cycle end if data=tdiurn(:,:,6) #ifdef TRACERS_AEROSOLS_Koch case ("SO4") ! sulfate in L=1 do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=trm(i,j,1,n_SO4) #ifdef TRACERS_HETCHEM * +trm(i,j,1,n_SO4_d1) * +trm(i,j,1,n_SO4_d2) * +trm(i,j,1,n_SO4_d3) #endif end do end do #endif #ifdef TRACERS_COSMO case ("7BEW") data=Be7w_acc case ("7BED") data=Be7d_acc #endif case default goto 10 end select kunit=kunit+1 C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) C**** write out !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) cycle C**** diags on fixed pressure levels or velocity 10 select case (namedd(k)(1:1)) case ("Z","R","T","Q") ! heights, rel/spec humidity or temp C**** get pressure level do kp=1,kgz_max if (namedd(k)(2:5) .eq. PMNAME(kp)) then kunit=kunit+1 select case (namedd(k)(1:1)) case ("Z") ! geopotential heights data=z_inst(kp,:,:) case ("R") ! relative humidity (wrt water) data=rh_inst(kp,:,:) case ("Q") ! specific humidity do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=rh_inst(kp,i,j)*qsat(t_inst(kp,i,j),lhe * ,PMB(kp)) end do end do case ("T") ! temperature (C) data=t_inst(kp,:,:) end select C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) C**** write out !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) cycle end if end do if (namedd(k)(2:4) .eq. "ALL") then do kp=1,kgz_max kunit=kunit+1 select case (namedd(k)(1:1)) case ("Z") ! geopotential heights data=z_inst(kp,:,:) case ("R") ! relative humidity (wrt water) data=rh_inst(kp,:,:) case ("Q") ! specific humidity do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=rh_inst(kp,i,j)*qsat(t_inst(kp,i,j),lhe * ,PMB(kp)) end do end do case ("T") ! temperature (C) data=t_inst(kp,:,:) end select C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) C**** write out !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do cycle end if case ("U","V","W","C") ! velocity/clouds on levels if (namedd(k)(2:4) .eq. "ALL") then kunit=kunit+1 do kp=1,LmaxSUBDD select case (namedd(k)(1:1)) case ("U") ! E-W velocity data=u(:,:,kp) case ("V") ! N-S velocity data=v(:,:,kp) case ("W") ! vertical velocity data=wsave(:,:,kp) C**** fix polar values for W only (calculated on tracer points) IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) case ("C") ! estimate of cloud optical depth do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=(1.-fss(kp,i,j))*taumc(kp,i,j)+fss(kp,i,j) * *tauss(kp,i,j) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) end select C**** write out !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do cycle end if C**** get model level do l=1,lm if (l.lt.10) then write(namel,'(I1)') l else write(namel,'(I2)') l end if if (trim(namedd(k)(2:5)) .eq. trim(namel)) then kunit=kunit+1 select case (namedd(k)(1:1)) case ("U") ! U velocity data=u(:,:,l) case ("V") ! V velocity data=v(:,:,l) case ("W") ! W velocity data=wsave(:,:,l) C**** fix polar values for W only (calculated on tracer points) IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) end select !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) cycle end if end do #ifdef TRACERS_SPECIAL_Shindell case ("O") ! Ox ozone tracer (ppmv) if (namedd(k)(3:5) .eq. "ALL") then kunit=kunit+1 do kp=1,LmaxSUBDD do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=1.d6*trm(i,j,kp,n_Ox)*mass2vol(n_Ox)/ * (am(kp,i,j)*dxyp(j)) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) C**** write out !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do cycle end if C**** get model level do l=1,lm if (l.lt.10) then write(namel,'(I1)') l else write(namel,'(I2)') l end if if (trim(namedd(k)(3:6)) .eq. trim(namel)) then kunit=kunit+1 do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=1d6*trm(i,j,l,n_Ox)*mass2vol(n_Ox)/ * (am(l,i,j)*dxyp(j)) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) C**** write out !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) cycle end if end do #endif #ifdef TRACERS_COSMO case ("B") ! Be7 tracer if (namedd(k)(2:4) .eq. "ALL") then kunit=kunit+1 do kp=1,LmaxSUBDD do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=1.d6*trm(i,j,kp,n_Be7)* mass2vol(n_Be7)/ * (am(kp,i,j)*dxyp(j)) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) C**** write out !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do cycle end if C**** get model level do l=1,lm if (l.lt.10) then write(namel,'(I1)') l else write(namel,'(I2)') l end if if (trim(namedd(k)(2:5)) .eq. trim(namel)) then kunit=kunit+1 do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=1d6*trm(i,j,l,n_Be7)*mass2vol(n_Be7)/ * (am(l,i,j)*dxyp(j)) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) C**** write out !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) cycle end if end do #endif #ifdef TRACERS_SPECIAL_O18 case ("D") ! HDO tracer (permil) if (namedd(k)(2:4) .eq. "ALL") then kunit=kunit+1 do kp=1,LmaxSUBDD do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=1d3*(trm(i,j,kp,n_HDO)/(trm(i,j,kp,n_water) * *trw0(n_HDO))-1.) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) C**** write out !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do cycle end if C**** get model level do l=1,lm if (l.lt.10) then write(namel,'(I1)') l else write(namel,'(I2)') l end if if (trim(namedd(k)(2:5)) .eq. trim(namel)) then kunit=kunit+1 do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=1d3*(trm(i,j,kp,n_HDO)/(trm(i,j,kp,n_water) * *trw0(n_HDO))-1.) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) C**** write out !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) cycle end if end do #endif end select C**** Additional diags select case (namedd(k)) #ifdef TRACERS_HETCHEM case ("SO2") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=trm(i,j,l,n_SO2) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("SO4") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=trm(i,j,l,n_SO4) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("SO4_d1") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=trm(i,j,l,n_SO4_d1) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("SO4_d2") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)= trm(i,j,l,n_SO4_d2) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("SO4_d3") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=trm(i,j,l,n_SO4_d3) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("Clay") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=trm(i,j,l,n_Clay) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("Silt1") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=trm(i,j,l,n_Silt1) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("Silt2") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=trm(i,j,l,n_Silt2) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("Silt3") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=trm(i,j,l,n_Silt3) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do #endif #ifdef CLD_AER_CDNC !for 3 hrly diagnostics case ("CTEM") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=ctem(l,i,j) ! cld temp (K) at cld top c write(6,*)"CTEM1",ctem(l,i,j),i,j,l,k end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("CL3D") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=cl3d(l,i,j) ! cld LWC (kg m-3) c write(6,*)"CL3D",cl3d(l,i,j),i,j,l,k end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("CI3D") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=ci3d(l,i,j) ! cld IWC (kg m-3) c write(6,*)"CI3D",ci3d(l,i,j),i,j,l,k end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("CD3D") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=cd3d(l,i,j) ! cld thickness (m) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("CLDSS") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=100.d0*cldss(l,i,j) ! Cld cover LS(%) c if(cldss(l,i,j).gt.0.d0) c & write(6,*)"CLDSS",cldss(l,i,j),i,j,l,k end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("CLDMC") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=100.d0*cldmc(l,i,j) ! Cld cover MC(%) end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("CDN3D") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=cdn3d(l,i,j) ! cld CDNC (cm^-3) c if(cdn3d(l,i,j).gt.20.d0)write(6,*)"CDN",cdn3d(l,i,j),i,j,l,k end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("CRE3D") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=1.d-6*cre3d(l,i,j) ! cld Reff (m) c if(cre3d(l,i,j).gt.2.d0)write(6,*)"CRe",cre3d(l,i,j),i,j,l,k end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("TAUSS") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=tauss(l,i,j) ! LS cld tau end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("TAUMC") kunit=kunit+1 do l=1,lm do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=taumc(l,i,j) ! MC cld tau end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end do case ("CLWP") !LWP (kg m-2) kunit=kunit+1 do j=J_0,J_1 do i=1,imaxj(j) data(i,j)=clwp(i,j) c if(clwp(i,j).ne.0.) write(6,*)"CLWP",data(i,j),i,j,k,kunit end do end do C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) C**** write out !call writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) #endif #if (defined TRACERS_DUST) || (defined TRACERS_MINERALS) ||\ (defined TRACERS_QUARZHEM) CASE ('DUEMIS') ! Dust emission flux [kg/m^2/s] kunit=kunit+1 DO n=1,Ntm_dust data(:,:)=dust_flux_glob(:,:,n) C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !CALL writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) END DO #ifdef TRACERS_DUST CASE ('DUEMIS2') ! Dust emission flux 2 (diag. var. only) [kg/m^2/s] kunit=kunit+1 DO n=1,Ntm_dust data(:,:)=dust_flux2_glob(:,:,n) C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !CALL writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) END DO #endif #ifdef TRACERS_DRYDEP CASE ('DUDEPTURB') ! Turb. deposition flux of dust tracers [kg/m^2/s] kunit=kunit+1 DO n=1,Ntm_dust n1=n_fidx+n-1 IF (dodrydep(n1)) THEN data(:,:)=depo_turb_glob(:,:,1,n1)+ & depo_turb_glob(:,:,2,n1)+depo_turb_glob(:,:,3,n1) & +depo_turb_glob(:,:,4,n1) C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !CALL writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) END IF END DO CASE ('DUDEPGRAV') ! Gravit. settling flux of dust tracers [kg/m^2/s] kunit=kunit+1 DO n=1,Ntm_dust n1=n_fidx+n-1 IF (dodrydep(n1)) THEN data(:,:)=depo_grav_glob(:,:,1,n1) & +depo_grav_glob(:,:,2,n1)+depo_grav_glob(:,:,3,n1) & +depo_grav_glob(:,:,4,n1) C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !CALL writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) END IF END DO #endif CASE ('DUDEPWET') ! Wet deposition flux of dust tracers [kg/m^2/s] kunit=kunit+1 DO n=1,Ntm_dust #ifdef TRACERS_WATER n1=n_fidx+n-1 IF (dowetdep(n1)) THEN DO j=J_0,J_1 data(:,j)=trprec(n1,:,j)*bydxyp(j)/Dtsrc END DO #else DO j=J_0,J_1 data(:,j)=trprec_dust(n,:,j)*bydxyp(j)/Dtsrc END DO #endif C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !CALL writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) #ifdef TRACERS_WATER END IF #endif END DO CASE ('DUTRS') ! Mixing ratio of dust tracers at surface [kg/kg] kunit=kunit+1 DO n=1,Ntm_dust n1=n_fidx+n-1 data(:,:)=trs_glob(:,:,1,n1)+trs_glob(:,:,2,n1) & +trs_glob(:,:,3,n1)+trs_glob(:,:,4,n1) C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !CALL writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) END DO CASE ('DULOAD') ! Dust load [kg/m^2] kunit=kunit+1 DO n=1,Ntm_dust n1=n_fidx+n-1 data(:,:)=0.D0 DO j=J_0,J_1 DO l=1,LmaxSUBDD data(:,j)=data(:,j)+trm(:,j,l,n1) END DO data(:,j)=data(:,j)*bydxyp(j) END DO C**** fix polar values IF(HAVE_SOUTH_POLE) data(2:im,1) =data(1,1) IF(HAVE_NORTH_POLE) data(2:im,jm)=data(1,jm) !CALL writei(iu_subdd(kunit),itime,data,im*jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) END DO #endif end select end do ! end of k=1,kdd loop c**** return end subroutine get_subdd subroutine write_data(data,kunit),4 !@sum write out subdd data array USE DOMAIN_DECOMP, only : GRID,GET,WRITEI_PARALLEL real*4 data(im,jm) integer kunit logical :: have_south_pole,have_north_pole call get(grid,have_south_pole=have_south_pole, & have_north_pole=have_north_pole) c**** fix polar values if(have_south_pole) data(2:im,1) =data(1,1) if(have_north_pole) data(2:im,jm)=data(1,jm) call writei_parallel(grid, & iu_subdd(kunit),nameunit(iu_subdd(kunit)),data,itime) end subroutine write_data end module subdaily subroutine ahourly 1,10 !@sum ahourly saves instantaneous variables at sub-daily frequency !@+ for diurnal cycle diagnostics !@auth Reha Cakmur/Jan Perlwitz USE MODEL_COM, only : u,v,t,p,q,jdate,jhour,ptop,sig USE CONSTANT, only : bygrav USE domain_decomp,ONLY : am_i_root,get,globalsum,grid USE GEOM, only : imaxj,dxyp,bydxyp USE DYNAMICS, only : phi,wsave,pek,byam USE rad_com,ONLY : cosz1,srnflb_save,trnflb_save,ttausv_save, & ttausv_cs_save USE diag_com,ONLY : adiurn_dust,ndiupt,ndiuvar,ijdd,adiurn #ifndef NO_HDIURN & ,hdiurn #endif #ifdef TRACERS_DUST * ,idd_u1,idd_v1,idd_uv1,idd_t1,idd_qq1,idd_p1,idd_w1,idd_phi1 * ,idd_sr1,idd_tr1,idd_load1,idd_conc1,idd_tau1,idd_tau_cs1 #endif #ifdef TRACERS_ON USE TRACER_COM, only : trm #ifdef TRACERS_DUST & ,Ntm_dust,n_clay #endif #endif IMPLICIT NONE INTEGER :: i,j,ih,ihm,kr,n,n1,nx REAL*8 :: psk INTEGER,PARAMETER :: lmax_dd2=11, n_idxd=14*lmax_dd2 INTEGER :: idxd(n_idxd) REAL*8 :: tmp(NDIUVAR) REAL*8,DIMENSION(n_idxd,grid%J_STRT_HALO:grid%J_STOP_HALO, & NDIUPT) :: diurn_partd REAL*8,DIMENSION(n_idxd, NDIUPT) :: diurnsumd C**** define local grid INTEGER J_0, J_1 C**** C**** Extract useful local domain parameters from "grid" C**** CALL get(grid, J_STRT=J_0, J_STOP=J_1) #ifdef TRACERS_DUST IF (adiurn_dust == 1) THEN diurn_partd=0.D0 idxd=(/ * (idd_u1+i-1,i=1,lmax_dd2), (idd_v1+i-1,i=1,lmax_dd2), * (idd_uv1+i-1,i=1,lmax_dd2), (idd_t1+i-1,i=1,lmax_dd2), * (idd_qq1+i-1,i=1,lmax_dd2), (idd_p1+i-1,i=1,lmax_dd2), * (idd_w1+i-1,i=1,lmax_dd2), (idd_phi1+i-1,i=1,lmax_dd2), * (idd_sr1+i-1,i=1,lmax_dd2), (idd_tr1+i-1,i=1,lmax_dd2), * (idd_load1+i-1,i=1,lmax_dd2), (idd_conc1+i-1,i=1,lmax_dd2), * (idd_tau1+i-1,i=1,lmax_dd2), (idd_tau_cs1+i-1,i=1,lmax_dd2) * /) END IF #endif ih=jhour+1 ihm=ih+(jdate-1)*24 !$OMP PARALLEL DO PRIVATE(i,j,kr,n,psk,n1,tmp) !$OMP* SCHEDULE(DYNAMIC,2) do j=j_0,j_1 do i=1,imaxj(j) psk=pek(1,i,j) do kr=1,ndiupt if(i.eq.ijdd(1,kr).and.j.eq.ijdd(2,kr)) then #ifdef TRACERS_DUST IF (adiurn_dust == 1) THEN tmp=0.D0 tmp(idd_u1:idd_u1+lmax_dd2-1)=u(i,j,1:lmax_dd2) tmp(idd_v1:idd_v1+lmax_dd2-1)=v(i,j,1:lmax_dd2) tmp(idd_uv1:idd_uv1+lmax_dd2-1)=sqrt(u(i,j,1:lmax_dd2)*u(i * ,j,1:lmax_dd2)+v(i,j,1:lmax_dd2)*v(i,j,1:lmax_dd2)) tmp(idd_t1:idd_t1+lmax_dd2-1)=t(i,j,1:lmax_dd2)*psk tmp(idd_qq1:idd_qq1+lmax_dd2-1)=q(i,j,1:lmax_dd2) tmp(idd_p1:idd_p1+lmax_dd2-1)=p(i,j)*sig(1:lmax_dd2)+ptop tmp(idd_w1:idd_w1+lmax_dd2-1)=wsave(i,j,1:lmax_dd2) tmp(idd_phi1:idd_phi1+lmax_dd2-1)=phi(i,j,1:lmax_dd2)*bygrav tmp(idd_sr1:idd_sr1+lmax_dd2-1)=srnflb_save(i,j,1:lmax_dd2) * *cosz1(i,j) tmp(idd_tr1:idd_tr1+lmax_dd2-1)=trnflb_save(i,j,1:lmax_dd2) DO n=1,Ntm_dust n1=n_clay+n-1 tmp(idd_load1:idd_load1+lmax_dd2-1) * =tmp(idd_load1:idd_load1+lmax_dd2-1)+trm(i,j * ,1:lmax_dd2,n1)*bydxyp(j) tmp(idd_conc1:idd_conc1+lmax_dd2-1) * =tmp(idd_conc1:idd_conc1+lmax_dd2-1)+trm(i,j * ,1:lmax_dd2,n1)*byam(1,i,j)*bydxyp(j) tmp(idd_tau1:idd_tau1+lmax_dd2-1)=tmp(idd_tau1:idd_tau1 * +lmax_dd2-1)+ttausv_save(i,j,n1,1:lmax_dd2) tmp(idd_tau_cs1:idd_tau_cs1+lmax_dd2-1) * =tmp(idd_tau_cs1:idd_tau_cs1+lmax_dd2-1) * +ttausv_cs_save(i,j,n1,1:lmax_dd2) END DO DIURN_partd(:,J,kr)=DIURN_partd(:,J,kr)+tmp(idxd(:)) END IF #endif endif enddo enddo enddo !$OMP END PARALLEL DO #ifdef TRACERS_DUST IF (adiurn_dust == 1) THEN CALL globalsum(grid,diurn_partd,diurnsumd) IF (AM_I_ROOT()) THEN adiurn(ih,idxd,:)=adiurn(ih,idxd,:)+diurnsumd #ifndef NO_HDIURN hdiurn(ihm,idxd,:)=hdiurn(ihm,idxd,:)+diurnsumd #endif END IF END IF #endif return end subroutine ahourly SUBROUTINE init_DIAG(istart,num_acc_files) 2,54 !@sum init_DIAG initializes the diagnostics !@auth Gavin Schmidt !@ver 1.0 USE CONSTANT, only : sday,kapa,undef USE MODEL_COM, only : lm,Itime,ItimeI,Itime0,pmtop,nfiltr,jhour * ,jdate,jmon,amon,jyear,jhour0,jdate0,jmon0,amon0,jyear0,idacc * ,ioread_single,xlabel,iowrite_single,iyear1,nday,dtsrc,dt * ,nmonav,ItimeE,lrunid,focean,pednl00,pmidl00,lm_req USE GEOM, only : imaxj USE SEAICE_COM, only : rsi USE LAKES_COM, only : flake USE DIAG_COM, only : TSFREZ => TSFREZ_loc USE DIAG_COM, only : NPTS, NAMDD, NDIUPT, IJDD, ISCCP_DIAGS USE DIAG_COM, only : monacc, acc_period, keyct, KEYNR, PLE USE DIAG_COM, only : PLM, p1000k, icon_AM, NOFM USE DIAG_COM, only : PLE_DN, icon_KE, NSUM_CON, IA_CON, SCALE_CON USE DIAG_COM, only : TITLE_CON, PSPEC, LSTR, NSPHER, KLAYER USE DIAG_COM, only : ISTRAT, kgz, pmb, kgz_max USE DIAG_COM, only : TF_DAY1, TF_LAST, TF_LKON, TF_LKOFF USE DIAG_COM, only : name_consrv, units_consrv, lname_consrv USE DIAG_COM, only : CONPT0, icon_MS, icon_TPE, icon_WM, icon_EWM USE diag_com,ONLY : adiurn_dust USE DIAG_LOC USE PARAM USE FILEMANAGER USE DOMAIN_DECOMP, only: GRID,GET,WRITE_PARALLEL IMPLICIT NONE integer, intent(in) :: istart,num_acc_files INTEGER I,J,L,K,KL,ioerr,months,years,mswitch,ldate,iu_AIC * ,jday0,jday,moff,kb,iu_ACC,l850,l300,l50 REAL*8 PLE_tmp CHARACTER CONPT(NPTS)*10 LOGICAL :: QCON(NPTS), T=.TRUE. , F=.FALSE. INTEGER :: J_0,J_1 !@var out_line local variable to hold mixed-type output for parallel I/O character(len=300) :: out_line CALL GET(GRID,J_STRT=J_0,J_STOP=J_1) C**** Initialise resolution dependent diagnostic parameters call diag_res call sync_param( "NAMDD", NAMDD, NDIUPT ) call sync_param( "IJDD", IJDD(1:2,1), 2*NDIUPT ) call sync_param( "isccp_diags",isccp_diags) call sync_param( "adiurn_dust",adiurn_dust) IF(ISTART.LT.1) THEN ! initialize for post-processing call getdte(Itime0,Nday,Iyear1,Jyear0,Jmon0,Jday0,Jdate0,Jhour0 * ,amon0) call getdte(Itime,Nday,Iyear1,Jyear,Jmon,Jday,Jdate,Jhour * ,amon) months=1 ; years=monacc(jmon0) ; mswitch=0 ; moff=0 ; kb=jmon0 do kl=jmon0+1,jmon0+11 k = kl if (k.gt.12) k=k-12 if (monacc(k).eq.years) then months=months+1 else if (monacc(k).ne.0) then C**** write(6,*) 'uneven period:',monacc CALL WRITE_PARALLEL(monacc, UNIT=6, format= & "('uneven period:',12I5)") call stop_model( 'uneven period', 255 ) end if if(monacc(k).ne.monacc(kb)) mswitch = mswitch+1 if(mswitch.eq.2) moff = moff+1 kb = k end do if (mswitch.gt.2) then C**** write(6,*) 'non-consecutive period:',monacc CALL WRITE_PARALLEL(monacc, UNIT=6, format= & "('non-consecutive period:',12I5)") call stop_model( 'non-consecutive period', 255 ) end if call aPERIOD (JMON0,JYEAR0,months,years,moff, acc_period,Ldate) if (num_acc_files.gt.1) then ! save the summed acc-file write(out_line,*) num_acc_files,' files are summed up' CALL WRITE_PARALLEL(TRIM(out_line), UNIT=6) keyct=1 ; KEYNR=0 XLABEL(128:132)=' ' XLABEL(120:132)=acc_period(1:3)//' '//acc_period(4:Ldate) C**** write(6,*) XLABEL CALL WRITE_PARALLEL(XLABEL, UNIT=6) call openunit(acc_period(1:Ldate)//'.acc'//XLABEL(1:LRUNID) * ,iu_ACC,.true.,.false.) call io_rsf (iu_ACC,Itime,iowrite_single,ioerr) call closeunit(iu_ACC) end if ItimeE = -1 close (6) open(6,file=acc_period(1:Ldate)//'.'//XLABEL(1:LRUNID)//'.PRT', * FORM='FORMATTED') END IF C**** Initialize certain arrays used by more than one print routine DO L=1,LM PLE(L) =pednl00(l+1) PLE_DN(L)=pednl00(l) PLM(L) =pmidl00(l) END DO PLM(LM+1:LM+LM_REQ)=pmidl00(lm+1:lm+lm_req) p1000k=1000.0**kapa C**** Initialise some local constants (replaces IFIRST constructions) C**** From DIAGA: DO L=1,LM LUPA(L)=L+1 LDNA(L)=L-1 END DO LDNA(1)=1 LUPA(LM)=LM C**** From DIAGB (PM, PMO are fixed, PL,PLO will vary) PM(1)=1200. ! ensures below surface for extrapolation DO L=2,LM+1 PL(L)=pednl00(l) PM(L)=pednl00(l) END DO DO L=1,LM PLO(L)=pmidl00(l) PMO(L)=.5*(PM(L)+PM(L+1)) END DO C**** From DIAG7A L850=LM L300=LM L50=LM DO L=LM-1,1,-1 PLE_tmp=.25*(PEDNL00(L)+2.*PEDNL00(L+1)+PEDNL00(L+2)) IF (PLE_tmp.LT.850.) L850=L IF (PLE_tmp.LT.300.) L300=L IF (PLE_tmp.LT.250.) JET=L IF (PLE_tmp.LT.50.) L50=L END DO C WRITE (6,888) JET CALL WRITE_PARALLEL(JET, UNIT=6, format= & "(' JET WIND LEVEL FOR DIAG',I3)") C 888 FORMAT (' JET WIND LEVEL FOR DIAG',I3) C**** WRITE (6,889) L850,L300,L50 WRITE (out_line,889) L850,L300,L50 CALL WRITE_PARALLEL(trim(out_line), UNIT=6) 889 FORMAT (' LEVELS FOR WIND WAVE POWER DIAG L850=',I3, * ' L300=',I3,' L50=',I3) LDEX(1)=L850 LDEX(2)=L300 LDEX(3)=L50 C**** Initialize conservation diagnostics C**** NCON=1:25 are special cases: Angular momentum and kinetic energy icon_AM=1 NOFM(:,icon_AM) = (/ 1, 8, 0, 0, 0, 0, 9,10, 0,11, 0, 0/) icon_KE=2 NOFM(:,icon_KE) = (/ 13,20,21, 0, 0, 0,22,23, 0,24, 0, 0/) NSUM_CON(1:25) = (/-1,-1,-1,-1,-1,-1,-1,12,12,12,12, 0, * -1,-1,-1,-1,-1,-1,-1,25,25,25,25,25, 0/) IA_CON(1:25) = (/12, 1, 1, 1, 1, 1, 1, 7, 8,10, 9,12, * 12, 1, 1, 1, 1, 1, 1, 7, 8, 8,10, 9,12/) SCALE_CON(1) = 1d-9 SCALE_CON((/2,3,4,5,6,7,8,9/))= 1d-2/DTSRC SCALE_CON(10) = 1d-2/(NFILTR*DTSRC) SCALE_CON(11) = 2d-2/SDAY SCALE_CON((/12,25/)) = 1. SCALE_CON(13) = 1d-3 SCALE_CON((/14,15,16,17,18,19,20,21,22/)) = 1d3/DTSRC SCALE_CON(23) = 1d3/(NFILTR*DTSRC) SCALE_CON(24) = 2d3/SDAY TITLE_CON(1:25) = (/ * ' INSTANTANE AM (10**9 J*S/M^2) ', * ' DELTA AM BY ADVECTION ', * ' DELTA AM BY CORIOLIS FORCE ', * ' DELTA AM BY PRESSURE GRAD ', * ' DELTA AM BY STRATOS DRAG ', * ' DELTA AM BY UV FILTER ', * ' DELTA AM BY GW DRAG ', * ' CHANGE OF AM BY DYNAMICS ', * ' CHANGE OF AM BY SURF FRIC+TURB ', * ' CHANGE OF AM BY FILTER ', * ' CHANGE OF AM BY DAILY RESTOR ', * ' SUM OF CHANGES (10**2 J/M^2) ', * '0INSTANTANEOUS KE (10**3 J/M^2) ', * ' DELTA KE BY ADVECTION ', * ' DELTA KE BY CORIOLIS FORCE ', * ' DELTA KE BY PRESSURE GRAD ', * ' DELTA KE BY STRATOS DRAG ', * ' DELTA KE BY UV FILTER ', * ' DELTA KE BY GW DRAG ', * ' CHANGE OF KE BY DYNAMICS ', * ' CHANGE OF KE BY MOIST CONVEC ', * ' CHANGE OF KE BY SURF + DC/TURB ', * ' CHANGE OF KE BY FILTER ', * ' CHANGE OF KE BY DAILY RESTOR ', * ' SUM OF CHANGES (10**-3 W/M^2) '/) name_consrv(1:25) = (/ * 'inst_AM ','del_AM_ADV','del_AM_COR','del_AM_PRE', * 'del_AM_STR','del_AM_UVF','del_AM_GWD','chg_AM_DYN' * ,'chg_AM_SUR','chg_AM_FIL','chg_AM_DAI','sum_chg_AM' * ,'inst_KE ','del_KE_ADV','del_KE_COR','del_KE_PRE' * ,'del_KE_STR','del_KE_UVF','del_KE_GWD','chg_KE_DYN' * ,'chg_KE_MOI','chg_KE_SUR','del_KE_FIL','chg_KE_DAI' * ,'sum_chg_KE'/) units_consrv(1) ="10**9 J*S/M^2" units_consrv(2:12) ="10**2 J/M^2" units_consrv(13) ="10**3 J/M^2" units_consrv(14:24)="10**-3 W/M^2" lname_consrv(1:25)=TITLE_CON(1:25) C**** To add a new conservation diagnostic: C**** i) Add 1 to NQUANT, and increase KCON in DIAG_COM.f C**** ii) Set up a QCON, and call SET_CON to allocate array numbers, C**** set up scales, titles, etc. The icon_XX index must be C**** declared in DIAG_COM.f for the time being C**** QCON denotes when the conservation diags should be done C**** 1:NPTS ==> DYN, COND, RAD, PREC, LAND, SURF, C**** FILTER,STRDG/OCEAN, DAILY, OCEAN1, OCEAN2, C**** iii) Write a conserv_XYZ routine that returns the zonal average C**** of your quantity C**** iv) Add a line to DIAGCA that calls conserv_DIAG (declared C**** as external) C**** v) Note that the conserv_XYZ routine, and call to SET_CON C**** should be in the driver module for the relevant physics C**** Set up atmospheric component conservation diagnostics CONPT=CONPT0 C**** Atmospheric mass QCON=(/ T, F, F, F, F, F, T, F, T, F, F/) CALL SET_CON(QCON,CONPT,"MASS ","(KG/M^2) ", * "(10**-8 KG/SM^2)",1d0,1d8,icon_MS) C**** Atmospheric total potential energy CONPT(8)="SURF+TURB" ; CONPT(6)="KE DISSIP" QCON=(/ T, T, T, F, F, T, T, T, F, F, F/) CALL SET_CON(QCON,CONPT,"TPE ","(10**5 J/M^2) ", * "(10**-2 W/M^2) ",1d-5,1d2,icon_TPE) C**** Atmospheric water mass CONPT(6)="SURF+TURB" QCON=(/ T, T, F, F, F, T, F, F, F, F, F/) CALL SET_CON(QCON,CONPT,"ATM WAT ","(10**-2 KG/M^2) ", * "(10**-8 KG/SM^2)",1d2,1d8,icon_WM) C**** Atmospheric water latent heat (at some point should include C**** sensible + potential energy associated with water mass as well) QCON=(/ T, T, F, F, F, T, F, F, F, F, F/) CALL SET_CON(QCON,CONPT,"ENRG WAT","(10**3 J/M^2) ", * "(10**-2 W/M^2) ",1d-3,1d2,icon_EWM) C**** Initialize layering for spectral diagnostics C**** add in epsilon=1d-5 to avoid roundoff mistakes KL=1 DO L=1,LM IF (PEDNL00(L+1)+1d-5.lt.PSPEC(KL) .and. * PEDNL00(L) +1d-5.gt.PSPEC(KL)) THEN IF (KL.eq.2) LSTR = L ! approx. 10mb height KL=KL+1 END IF KLAYER(L)=4*(KL-1)+1 END DO IF (KL*4 .gt. NSPHER) THEN C**** WRITE(6,*) "Inconsistent definitions of stratosphere:" CALL WRITE_PARALLEL("Inconsistent definitions of stratosphere:" & ,UNIT=6) C**** WRITE(6,*) "Adjust PSPEC, ISTRAT so that KL*4 = NSPHER" CALL WRITE_PARALLEL("Adjust PSPEC, ISTRAT so that KL*4 = NSPHER" & ,UNIT=6) C**** WRITE(6,*) "ISTRAT,PSPEC,NSPHER,KL=",ISTRAT,PSPEC,NSPHER,KL WRITE(out_line,*) "ISTRAT,PSPEC,NSPHER,KL=", & ISTRAT,PSPEC,NSPHER,KL CALL WRITE_PARALLEL(trim(out_line), UNIT=6) call stop_model( * "Stratospheric definition problem for spectral diags.",255) END IF C**** Calculate the max number of geopotential heights do k=1,kgz if (pmb(k).le.pmtop) exit kgz_max = k end do C**** write(6,'(a)') " Geopotential height diagnostics at (mb): " CALL WRITE_PARALLEL(" Geopotential height diagnostics at (mb): ", & UNIT=6) C***** write(6,'(20F9.3)') PMB(1:kgz_max) CALL WRITE_PARALLEL(PMB(1:kgz_max), UNIT=6, format="(20F9.3)") c**** Initialize acc-array names, units, idacc-indices call def_acc C**** Ensure that diagnostics are reset at the beginning of the run IF (Itime.le.ItimeI .and. ISTART.gt.0) THEN call getdte(Itime,Nday,Iyear1,Jyear,Jmon,Jday,Jdate,Jhour * ,amon) CALL reset_DIAG(0) C**** Initiallise ice freeze diagnostics at beginning of run DO J=J_0,J_1 DO I=1,IMAXJ(J) TSFREZ(I,J,TF_DAY1)=365. TSFREZ(I,J,TF_LAST)=365. IF (FOCEAN(I,J)+FLAKE(I,J).gt.0) then IF (RSI(I,J).gt.0) then TSFREZ(I,J,TF_LKON) = JDAY-1 TSFREZ(I,J,TF_LKOFF) = JDAY ELSE TSFREZ(I,J,TF_LKON) = JDAY TSFREZ(I,J,TF_LKOFF) = undef END IF ELSE TSFREZ(I,J,TF_LKON) = undef TSFREZ(I,J,TF_LKOFF) = undef END IF END DO END DO CALL daily_DIAG END IF RETURN END SUBROUTINE init_DIAG SUBROUTINE reset_DIAG(isum) 4,8 !@sum reset_DIAG resets/initializes diagnostics !@auth Original Development Team !@ver 1.0 USE MODEL_COM, only : Itime,iyear1,nday,kradia, * Itime0,jhour0,jdate0,jmon0,amon0,jyear0,idacc,u USE DIAG_COM USE PARAM #ifdef TRACERS_ON USE TRDIAG_COM, only: taijln, TAIJLN_loc, taijn, TAIJN_loc, * taijs, TAIJS_loc, TAJLN, TAJLN_loc, TAJLS, TAJLS_loc, * TCONSRV,TCONSRV_loc #endif USE DOMAIN_DECOMP, only: grid, CHECKSUM IMPLICIT NONE INTEGER :: isum !@var isum if =1 preparation to add up acc-files INTEGER jd0 IDACC(1:12)=0 if (kradia.gt.0) then AFLX_ST = 0. if (isum.eq.1) return go to 100 end if AJ_loc=0 ; AREGJ_loc=0 APJ_loc=0 ; AJL_loc=0 ; ASJL_loc=0 ; AIJ_loc=0 AIL=0 ; ENERGY=0 ; CONSRV_loc=0 SPECA=0 ; ATPE=0 ; WAVE=0 ; AJK_loc=0 ; AIJK_loc=0 #ifndef NO_HDIURN HDIURN=0 #endif ADIURN=0 ; AISCCP=0 #ifdef TRACERS_ON TAJLN=0. ; TAJLS=0. ; TCONSRV=0. TAJLN_loc=0. ; TAJLS_loc=0. ; TCONSRV_loc=0. TAIJLN=0. ; TAIJN=0. ; TAIJS=0. TAIJLN_loc=0. ; TAIJN_loc=0. ; TAIJS_loc=0. #endif call reset_ODIAG(isum) ! ocean diags if required call reset_icdiag ! ice dynamic diags if required if (isum.eq.1) then ! prepare to add up acc-files AJ=0 ; APJ=0 ; AJL=0 ; ASJL=0 ; AIJ=0 ; AJK=0 ; AIJK=0 return end if AIJ_loc(:,:,IJ_TMNMX)=1000. ; IDACC(12)=1 CALL EPFLXI (U) ! strat 100 Itime0=Itime call getdte(Itime0,Nday,Iyear1,Jyear0,Jmon0,Jd0, * Jdate0,Jhour0,amon0) RETURN END SUBROUTINE reset_DIAG SUBROUTINE daily_DIAG 2,9 !@sum daily_DIAG resets diagnostics at beginning of each day !@auth Original Development Team !@ver 1.0 USE CONSTANT, only : undef USE MODEL_COM, only : im,jm,jday,focean USE GEOM, only : imaxj USE SEAICE_COM, only : rsi USE LAKES_COM, only : flake USE GHY_COM, only : fearth USE DIAG_COM, only : aij=>aij_loc * ,ij_lkon,ij_lkoff,ij_lkice,tsfrez=>tsfrez_loc,tdiurn * ,tf_lkon,tf_lkoff,tf_day1,tf_last USE DOMAIN_DECOMP, only : GRID,GET IMPLICIT NONE INTEGER I,J INTEGER :: J_0, J_1 CALL GET(GRID,J_STRT=J_0,J_STOP=J_1) C**** INITIALIZE SOME ARRAYS AT THE BEGINNING OF SPECIFIED DAYS IF (JDAY.EQ.32) THEN DO J=MAX(J_0,1+JM/2),MIN(J_1,JM) DO I=1,IM TSFREZ(I,J,TF_DAY1)=JDAY END DO END DO DO J=MAX(J_0,1),MIN(J_1,JM/2) DO I=1,IM TSFREZ(I,J,TF_LAST)=JDAY END DO END DO ELSEIF (JDAY.EQ.213) THEN DO J=MAX(J_0,1),MIN(J_1,JM/2) DO I=1,IM TSFREZ(I,J,TF_DAY1)=JDAY END DO END DO END IF C**** set and initiallise freezing diagnostics C**** Note that TSFREZ saves the last day of no-ice and some-ice. C**** The AIJ diagnostics are set once a year (zero otherwise) DO J=J_0,J_1 DO I=1,IMAXJ(J) IF (J.le.JM/2) THEN C**** initiallise/save South. Hemi. on Feb 28 IF (JDAY.eq.59 .and. TSFREZ(I,J,TF_LKOFF).ne.undef) THEN AIJ(I,J,IJ_LKICE)=1. AIJ(I,J,IJ_LKON) =MOD(NINT(TSFREZ(I,J,TF_LKON)) +307,365) AIJ(I,J,IJ_LKOFF)=MOD(NINT(TSFREZ(I,J,TF_LKOFF))+306,365) * +1 IF (RSI(I,J).gt.0) THEN TSFREZ(I,J,TF_LKON) = JDAY-1 ELSE TSFREZ(I,J,TF_LKOFF) = undef END IF END IF ELSE C**** initiallise/save North. Hemi. on Aug 31 C**** Note that for continuity across the new year, the julian days C**** are counted from Sep 1 (NH only). IF (JDAY.eq.243 .and. TSFREZ(I,J,TF_LKOFF).ne.undef) THEN AIJ(I,J,IJ_LKICE)=1. AIJ(I,J,IJ_LKON) =MOD(NINT(TSFREZ(I,J,TF_LKON)) +123,365) AIJ(I,J,IJ_LKOFF)=MOD(NINT(TSFREZ(I,J,TF_LKOFF))+122,365) * +1 IF (RSI(I,J).gt.0) THEN TSFREZ(I,J,TF_LKON) = JDAY-1 ELSE TSFREZ(I,J,TF_LKOFF) = undef END IF END IF END IF C**** set ice on/off days IF (FOCEAN(I,J)+FLAKE(I,J).gt.0) THEN IF (RSI(I,J).eq.0.and.TSFREZ(I,J,TF_LKOFF).eq.undef) * TSFREZ(I,J,TF_LKON)=JDAY IF (RSI(I,J).gt.0) TSFREZ(I,J,TF_LKOFF)=JDAY END IF END DO END DO C**** INITIALIZE SOME ARRAYS AT THE BEGINNING OF EACH DAY DO J=J_0,J_1 DO I=1,IM TDIURN(I,J,1)= 1000. TDIURN(I,J,2)=-1000. TDIURN(I,J,3)= 1000. TDIURN(I,J,4)=-1000. TDIURN(I,J,5)= 0. TDIURN(I,J,6)=-1000. TDIURN(I,J,7)=-1000. TDIURN(I,J,8)=-1000. TDIURN(I,J,9)= 1000. IF (FEARTH(I,J).LE.0.) THEN TSFREZ(I,J,TF_DAY1)=365. TSFREZ(I,J,TF_LAST)=365. END IF END DO END DO END SUBROUTINE daily_DIAG SUBROUTINE SET_CON(QCON,CONPT,NAME_CON,INST_UNIT,SUM_UNIT,INST_SC 23,8 * ,CHNG_SC,ICON) !@sum SET_CON assigns conservation diagnostic array indices !@auth Gavin Schmidt !@ver 1.0 USE CONSTANT, only : sday USE MODEL_COM, only : dtsrc,nfiltr USE DIAG_COM, only : kcon,nquant,npts,title_con,scale_con,nsum_con * ,nofm,ia_con,kcmx,ia_d5d,ia_d5s,ia_filt,ia_12hr,name_consrv * ,lname_consrv,units_consrv USE DOMAIN_DECOMP, only : WRITE_PARALLEL IMPLICIT NONE !@var QCON logical variable sets where conservation diags are saved LOGICAL, INTENT(IN),DIMENSION(NPTS) :: QCON !@var CONPT names for points where conservation diags are saved CHARACTER*10, INTENT(IN),DIMENSION(NPTS) :: CONPT !@var INST_SC scale for instantaneous value REAL*8, INTENT(IN) :: INST_SC !@var CHNG_SC scale for changes REAL*8, INTENT(IN) :: CHNG_SC !@var NAME_CON name of conservation quantity CHARACTER*8, INTENT(IN) :: NAME_CON !@var sname name of conservation quantity (no spaces) CHARACTER*8 :: sname !@var INST_UNIT string for unit for instant. values CHARACTER*16, INTENT(IN) :: INST_UNIT !@var SUM_UNIT string for unit for summed changes CHARACTER*16, INTENT(IN) :: SUM_UNIT !@var ICON index for the conserved quantity INTEGER, INTENT(OUT) :: ICON !@var out_line local variable to hold mixed-type output for parallel I/O character(len=300) :: out_line INTEGER NI,NM,NS,N,k INTEGER, SAVE :: NQ = 2 ! first 2 special cases AM + KE NQ=NQ+1 IF (NQ.gt.NQUANT) THEN C**** WRITE(6,*) "Number of conserved quantities larger than NQUANT" C**** * ,NQUANT,NQ WRITE(out_line,*) * "Number of conserved quantities larger than NQUANT" * ,NQUANT,NQ CALL WRITE_PARALLEL(trim(out_line), UNIT=6) call stop_model("Change NQUANT in diagnostic common block",255) END IF C**** remove spaces in NAME_CON for netcdf names sname=TRIM(NAME_CON) do k=1,len_trim(NAME_CON) if (sname(k:k).eq." ") sname(k:k)="_" end do C**** NI=KCMX+1 NOFM(1,NQ) = NI TITLE_CON(NI) = "0INSTANT "//TRIM(NAME_CON)//" "//TRIM(INST_UNIT) SCALE_CON(NI) = INST_SC name_consrv(NI) ="inst_"//sname lname_consrv(NI) = "INSTANT "//TRIM(NAME_CON) units_consrv(NI) = INST_UNIT IA_CON(NI) = 12 NM=NI DO N=1,NPTS IF (QCON(N)) THEN NM = NM + 1 NOFM(N+1,NQ) = NM TITLE_CON(NM) = " CHANGE OF "//TRIM(NAME_CON)//" BY "// * CONPT(N) name_consrv(NM) ="chg_"//trim(sname)//"_"//TRIM(CONPT(N)(1:3)) lname_consrv(NM) = TITLE_CON(NM) units_consrv(NM) = SUM_UNIT SELECT CASE (N) CASE (1) SCALE_CON(NM) = CHNG_SC/DTSRC IA_CON(NM) = ia_d5d CASE (2,3,4,5,6,8,10,11) SCALE_CON(NM) = CHNG_SC/DTSRC IA_CON(NM) = ia_d5s CASE (7) SCALE_CON(NM) = CHNG_SC/(NFILTR*DTSRC) IA_CON(NM) = ia_filt CASE (9) SCALE_CON(NM) = CHNG_SC*2./SDAY IA_CON(NM) = ia_12hr END SELECT ELSE NOFM(N+1,NQ) = 0 END IF END DO NS=NM+1 IF (NS.gt.KCON) THEN C**** WRITE(6,*) "KCON not large enough for extra conserv diags", C**** * KCON,NI,NM,NQ,NS,NAME_CON WRITE(out_line,*) * "KCON not large enough for extra conserv diags", * KCON,NI,NM,NQ,NS,NAME_CON CALL WRITE_PARALLEL(trim(out_line), UNIT=6) call stop_model("Change KCON in diagnostic common block",255) END IF TITLE_CON(NS) = " SUM OF CHANGES "//TRIM(SUM_UNIT) name_consrv(NS) ="sum_chg_"//trim(sname) lname_consrv(NS) = " SUM OF CHANGES OF "//TRIM(NAME_CON) units_consrv(NS) = SUM_UNIT SCALE_CON(NS) = 1. IA_CON(NS) = 12 NSUM_CON(NI) = -1 NSUM_CON(NI+1:NS-1) = NS NSUM_CON(NS) = 0 KCMX=NS ICON=NQ RETURN C**** END SUBROUTINE set_con SUBROUTINE UPDTYPE 5,7 !@sum UPDTYPE updates FTYPE array to ensure correct budget diagnostics !@auth Gavin Schmidt USE MODEL_COM, only : im,jm,focean,flice,itocean * ,itoice,itlandi,itearth,itlake,itlkice,ftype USE GEOM, only : imaxj USE SEAICE_COM, only : rsi USE LAKES_COM, only : flake USE GHY_COM, only : fearth USE DOMAIN_DECOMP, only : GRID,GET IMPLICIT NONE INTEGER I,J INTEGER :: J_0,J_1 CALL GET(GRID,J_STRT=J_0,J_STOP=J_1) DO J=J_0,J_1 DO I=1,IMAXJ(J) FTYPE(ITOICE ,I,J)=FOCEAN(I,J)*RSI(I,J) FTYPE(ITOCEAN,I,J)=FOCEAN(I,J)-FTYPE(ITOICE,I,J) FTYPE(ITLKICE,I,J)=FLAKE(I,J)*RSI(I,J) FTYPE(ITLAKE ,I,J)=FLAKE(I,J)-FTYPE(ITLKICE,I,J) C**** set land components of FTYPE array. Summation is necessary for C**** cases where Earth and Land Ice are lumped together FTYPE(ITLANDI,I,J)=0. FTYPE(ITEARTH,I,J)=FEARTH(I,J) FTYPE(ITLANDI,I,J)=FTYPE(ITLANDI,I,J)+FLICE(I,J) END DO END DO RETURN C**** END SUBROUTINE UPDTYPE