diff --git a/model/ATM_DRV.f b/model/ATM_DRV.f index b9e27bf..0cd3133 100644 --- a/model/ATM_DRV.f +++ b/model/ATM_DRV.f @@ -858,6 +858,9 @@ subroutine alloc_drv_atm() #ifdef BIOGENIC_EMISSIONS call alloc_biogenic_emis(grid) #endif +#if (defined CALCULATE_LIGHTNING) || (defined TRACERS_SPECIAL_Shindell) + call alloc_lightning(grid) +#endif #if (defined TRACERS_ON) || (defined TRACERS_OCEAN) || (defined TRACERS_WATER) call alloc_tracer_com(grid) #ifdef TRACERS_DRYDEP @@ -867,9 +870,6 @@ subroutine alloc_drv_atm() call alloc_tracer_special_lerner_com(grid) call alloc_linoz_chem_com(grid) #endif -#if (defined CALCULATE_LIGHTNING) || (defined TRACERS_SPECIAL_Shindell) - call alloc_lightning(grid) -#endif #ifdef TRACERS_SPECIAL_Shindell call alloc_trchem_shindell_com(grid) call alloc_tracer_sources(grid) @@ -1727,6 +1727,10 @@ subroutine accum_subdd_atm use atm_com, only: u,v,t,q,qcl,qci, pdsig,pmid,pedn,pk, & ualij,valij, zatmo,gz, wsave, ma,masum & ,ptropo,ltropo +#ifdef GCAP + use atm_com, only : MWs + use geom, only : byaxyp +#endif use domain_decomp_atm, only : grid,get=>getdomainbounds use resolution, only : lm,mtop USE GEOM, only: imaxj @@ -1845,6 +1849,29 @@ subroutine accum_subdd_atm sddarr2d(i,j) = t(i,j,ltropo(i,j))*pk(ltropo(i,j),i,j) enddo; enddo call inc_subdd(subdd,k,sddarr2d) + +#ifdef GCAP + case ('HFLUX') ! Based on shflux + sddarr2d = -atmsrf%sensht(:,:)/dtsrc ! Note: sign change for consistency with MERRA-2 + call inc_subdd(subdd,k,sddarr2d) + case ('EFLUX') + sddarr2d = -atmsrf%latht(:,:)/dtsrc ! Note: sign change for consistency with MERRA-2 + call inc_subdd(subdd,k,sddarr2d) + case ('SLP') + do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr2d(i,j) = + & slp(pedn(1,i,j),atmsrf%tsavg(i,j),bygrav*zatmo(i,j))*100. ! Convert to Pa + enddo; enddo + call inc_subdd(subdd,k,sddarr2d) + case ('PBLH') + call inc_subdd(subdd,k,atmsrf%dblavg) + case( 'TROPPT' ) + call inc_subdd(subdd,k,ptropo*100) ! Convert to Pa + case ('PS') + sddarr2d = pedn(1,:,:)*100. ! Convert to Pa + call inc_subdd(subdd,k,sddarr2d) +#endif + end select enddo enddo @@ -1951,6 +1978,34 @@ subroutine accum_subdd_atm sddarr(:,:,1:lm-1) = wsave sddarr(:,:,lm) = 0. call inc_subdd(subdd,k,sddarr) +#ifdef GCAP + case ('T') + do l=1,lmaxsubdd; do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr(i,j,l) = t(i,j,l)*pk(l,i,j) + enddo; enddo; enddo + call inc_subdd(subdd,k,sddarr) + case ('QV') + call inc_subdd(subdd,k,q) + case ('OMEGA') +! Downward pressure velocity +!**** omega(Pa/s) = MWs(mb*m^2) * byAXYP(1/m^2) * 100(Pa/mb) / DTSRC(s) + sddarr(:,:,:) = 0 + do l=1,lmaxsubdd; do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr(i,j,l) = MWs(I,J,L)*byaxyp(I,J)*100/DTSRC + enddo; enddo; enddo + call inc_subdd(subdd,k,sddarr) + case ('RH') + do l=1,lmaxsubdd; do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr(i,j,l) = + & q(i,j,l)/QSAT(t(i,j,l)*pk(l,i,j),LHE,pmid(l,i,j)) ! Leave as fraction + enddo; enddo; enddo + call inc_subdd(subdd,k,sddarr) + case ('U') + call inc_subdd(subdd,k,ualij,jdim=3) + case ('V') + call inc_subdd(subdd,k,valij,jdim=3) +#endif + end select enddo enddo diff --git a/model/CLOUDS2.F90 b/model/CLOUDS2.F90 index 081050b..993614e 100644 --- a/model/CLOUDS2.F90 +++ b/model/CLOUDS2.F90 @@ -24,6 +24,7 @@ module CLOUDS ,mb2kg use RESOLUTION, only : lm USE ATM_COM, only : pdsigl00 + USE GEOM, only : DXYP use MODEL_COM, only : dtsrc,itime ! ,coupled_chem use TimeConstants_mod, only: SECONDS_PER_HOUR use CLOUDS_COM, only : ncol @@ -238,6 +239,9 @@ module CLOUDS !@var CLDSSL large-scale cloud cover !@var SM,QM Vertical profiles of (T/p**kappa)*AIRM, q*AIRM real*8, dimension(LM) :: SSHR,DCTEI,TAUSSL,CLDSSL,TAUSSLIP +#ifdef GCAP + real*8, dimension(LM) :: TAUSSL3D,CLDSSL3D +#endif real*8, dimension(LM) :: SM,QM real*8, dimension(NMOM,LM) :: SMOM,QMOM,SMOMMC,QMOMMC, SMOMLS,QMOMLS @@ -403,6 +407,30 @@ module CLOUDS !@var ccp_cosp Mixing ratio of convective precipitation [kg/kg] REAL*8 ccp_cosp(LM+1) #endif +#ifdef GCAP + real*8 mc_up_mf(LM+1,2), & ! Plume updraft mass flux along edge (kg/m2/s) + mc_dd_mf(LM+1,2), & ! Plume downdraft mass flux alonge edge (kg/m2/s) + mc_up_ent(LM,2), & ! Plume updraft entrainment flux (kg/m2/s) + mc_up_det(LM,2), & ! Plume updraft detrainment flux (kg/m2/s) + mc_dd_ent(LM,2), & ! Plume downdraft entrainment flux (kg/m2/s) + mc_dd_det(LM,2), & ! Plume downdraft detrainment flux (kg/m2/s) + mc_dqcond(LM,2), & ! Moist convection condensation rate (kg H2O/kg dry air in cell/s) + mc_dqevap(LM,2), & ! Moist convection evaporation rate (kg H2O/kg dry air in cell/s) + ls_dqcond(LM), & ! Large-scale condensation rate (kg H2O/kg dry air in cell/s) + ls_dqevap(LM), & ! Large-scale evaporation rate (kg H2O/kg dry air in cell/s) + mc_pflx_i(LM+1), & ! Moist convection liquid precip flux along edge (kg/m2/s) + mc_pflx_l(LM+1), & ! Moist convection ice precip flux along edge (kg/m2/s) + ls_pflx_i(LM+1), & ! Large-scale liquid precip flux along edge (kg/m2/s) + ls_pflx_l(LM+1) ! Large-scale ice precip flux along edge (kg/m2/s) + integer mc_counts(LM,2,6), mc_plumes ! Counts to make sure we're not missing anything + ! 1 = updraft mass flux + ! 2 = updraft entrain + ! 3 = updraft detrain + ! 4 = downdraft mass flux + ! 5 = downdraft entrain + ! 6 = downdraft detrain + real*8 :: dZ(LM), dZm(LM+1) +#endif #endif contains @@ -581,7 +609,7 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) CDHM,CDHSUM,CLDM,CLDREF,CDHSUM1,CONTCE,CDHDRT,CONDMU, & ! DCW,DCG,DCI,DWCU,DMSE,DFP,DQSUM,DMMIX,DQ,DMSE1,DQSUM1, & - DDRAFT,DELTA,DQEVP,DDRUP,DDROLD, & + DDRAFT,DELTA,DQEVP,DDRUP,DDROLD, DP, & EPLUME,ETADN,ETAL1,EVPSUM,EDRAFT, & ! FLAMW,FLAMG,FLAMI,FG,FI,FMC1,FCLW, & @@ -604,6 +632,10 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) TADJ,TEMWM,TEM,TIG,TNX1,TP,TVP,TOLD,TOLD1,TTURB,TRATIO, & UMTEMP,VMTEMP,VT, & WMDN,WMUP,WMEDG,WMIX,W2TEM,WTEM,WCONST,WCUFRZ,WORK,WMAX,WV +#ifdef GCAP + real*8, dimension(LM) :: bykg + real*8 :: ccm_ref, kg_plume +#endif ! ! ! ********* SCALAR DEFINITIONS ********* @@ -886,7 +918,38 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) ccp_cosp(:) = 0. ccl_cosp(:) = 0. #endif +#ifdef GCAP + ! Reset arrays that track the convective mass fluxes + mc_up_mf(:,:) = 0d0 + mc_dd_mf(:,:) = 0d0 + mc_up_ent(:,:) = 0d0 + mc_up_det(:,:) = 0d0 + mc_dd_ent(:,:) = 0d0 + mc_dd_det(:,:) = 0d0 + mc_counts(:,:,:) = 0 + mc_plumes = 0 + ! Reset arrays that track the condensation/evaporation fluxes + mc_dqcond(:,:) = 0d0 + mc_dqevap(:,:) = 0d0 + mc_pflx_l(:) = 0d0 + mc_pflx_i(:) = 0d0 + ! Calculate layer thickness [m] with hypsometric equation + DO L=1,LM + dZ(L) = rgas * TL(L) * bygrav * log(ple(L)/ple(L+1)) + ENDDO + ! Calculate distance layer midpoints [m] (with special treatment at surface) for + ! precipitation flux calculation + dZm(1) = 0.5 * dZ(1) + DO L=2,LM+1 + dZm(L) = 0.5 * ( dZm(L-1) + dZm(L) ) + ENDDO + ! Calculate inverse mass in each layer (1/kg) + DO L=1,LM + bykg(L) = grav / ( 100. * ( ple(L) - ple(L+1) ) * DXYP(j_debug) ) + ENDDO #endif +#endif + !**** SAVE ORIG PROFILES SMOLD(:) = SM(:) SMOMOLD(:,:) = SMOM(:,:) @@ -1063,6 +1126,12 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) If (NPPL==2 .and. (.not.MC1 .or. MCCONT.lt.2)) Cycle AREA_PARTITION !**** Convection rose 2 or more layers for first AREA PARTITION +#ifdef GCAP + ! Reset aggregators + mc_dqevap(:,IC) = 0. + mc_dqcond(:,IC) = 0. +#endif + #ifdef TRACERS_ON If (NPPL==2) then call reset_tracer_work_arrays (lmin,lmax) @@ -1215,6 +1284,10 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) CLOUD_TOP: do L=LMIN+1,LM + ! Reset some scalars used by diagnostics + eplume = 0 + delta = 0 + !**** !**** TRIGGERING CONDITIONS FOR MOIST CONVECTION !**** @@ -1328,7 +1401,16 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) !**** CONDENSE VAPOR IN THE PLUME AND ADD LATENT HEAT call get_dq_cond(smp,qmp,plk(l),mplume,lhx,pl(l),dqsum,fqcond) - +#ifdef GCAP + ! Save rainwater production rate (kg/kg/s) + ! Convert dqsum ( kg H2O / kg plume * hPa plume ) to ( kg H2O / kg air / s ) + ! This may be overwritten in NPPL=2 + ! dqsum is in kg H2O/kg plume mass * hPa plume mass + if ( lhx .eq. lhe ) then + kg_plume = 1d2 * mplume * FMC1 * DXYP(j_debug) * bygrav + mc_dqcond(L,IC) = ( dqsum / mplume ) * kg_plume * bykg(L) * bydtsrc + endif +#endif if(DQSUM.gt.0. .and. QMP.gt.teeny) then QMOMP(xymoms) = QMOMP(xymoms)*(1.-FQCOND) SMP=SMP+SLH*DQSUM/PLK(L) @@ -1346,6 +1428,14 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) CDHEAT(L)=SLH*COND(L) ! calculate CDHEAT before add CONDV CDHSUM=CDHSUM+CDHEAT(L) COND(L)=COND(L)+CONDV(L-1) ! add in the vertical transported COND +#ifdef GCAP + ! Convert (kg/kg)*mb to kg/m2/s; flip sign since in updraft + IF ( TL(L-1) .gt. TF ) THEN + mc_pflx_l(L) = mc_pflx_l(L) - ( ( CONDV(L-1) / CCM(L-1) ) * ( PL(L) * 1d2 / ( RGAS * TL(L) ) ) ) / dZm(L) + ELSE + mc_pflx_i(L) = mc_pflx_i(L) - ( ( CONDV(L-1) / CCM(L-1) ) * ( PL(L) * 1d2 / ( RGAS * TL(L) ) ) ) / dZm(L) + ENDIF +#endif if (VLAT(L-1).ne.VLAT(L)) then SMP=SMP-(VLAT(L-1)-VLAT(L))*CONDV(L-1)*BYSHA/PLK(L) !phase change of uplifted condensate @@ -1708,6 +1798,15 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) DVM(K,L)=DVM(K,L)-V_0(K,L)*EPLUME-VMTEMP end do +#ifdef GCAP + !==================================================================== + ! Archive mass that entrains into the upward plume (kg/m2/s) + !==================================================================== + ! Done with entrainment in upward plume, so archive here. + ! Note: Do not accumulate here, as EPLUME may be overwritten in NPPL=2 + mc_up_ent(L,IC) = 100.*EPLUME*bygrav/dtsrc +#endif + end if ! big enough check end if ! non-zero entrainment check @@ -1730,9 +1829,17 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) DETALL(L,IC,LM) = DET(L)*100.*1000. endif #endif +#ifdef GCAP + !========================================================================= + ! Archive mass that detrains from upward plumes (occurs at cloud tops; kg/m2/s) + !========================================================================= + ! MPLUME is the mass of the plume in mb + ! DELTA is the fraction that detrains from the plume into the layer + mc_up_det(L,IC) = mc_up_det(L,IC) + 100.*DELTA*MPLUME*bygrav/dtsrc +#endif - DM(L)=DM(L)+DELTA*MPLUME - MPLUME=MPLUME*(1.D0-DELTA) + DM(L)=DM(L)+DELTA*MPLUME ! Mass that leaves plume + MPLUME=MPLUME*(1.D0-DELTA) ! Mass that remains in plume DSM(L)= DSM(L)+DELTA*SMP SMP = SMP *(1.-DELTA) DSMOM(xymoms,L)=DSMOM(xymoms,L)+DELTA*SMOMP(xymoms) @@ -1915,8 +2022,8 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) if(CONDP(L).gt.CONDP1(L)) CONDP(L)=CONDP1(L) CONDV(L)=COND(L)-CONDP1(L) ! part of COND transported up COND(L)=COND(L)-CONDV(L) ! CONDP1(L) - TAUMC1(L)=TAUMC1(L)-CONDV(L)*FMC1 - + TAUMC1(L)=TAUMC1(L)-CONDV(L)*FMC1 + #ifdef TRACERS_WATER FQCONDV=CONDV(L)/((COND(L)+CONDV(L))+teeny) TRCONDV(:,L)=FQCONDV*TRCOND(:,L) @@ -2015,10 +2122,15 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) !**** LOOP FROM TOP DOWN OVER POSSIBLE DOWNDRAFTS !**** DOWNDRAFT: do L=LDRAFT,1,-1 + LHX=VLAT(L) ! LHX consistency SLH=LHX*BYSHA TNX1=SMDN*PLK(L)/DDRAFT ! save for tracers + ! Reset some scalars used in diagnostics + ! Downdrafts happen outside of the NPPL loop. + eplume = 0. + call get_dq_evap(smdn,qmdn,plk(l),ddraft,lhx,pl(l),cond(l), dqsum,fqcond1) !**** EVAPORATE CONVECTIVE CONDENSATE IN DOWNDRAFT AND UPDATE DOWNDRAFT @@ -2029,6 +2141,13 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) if(DQEVP.gt.SMDN*PLK(L)/SLH) DQEVP=SMDN*PLK(L)/SLH if (L.lt.LMIN) DQEVP=0. +#ifdef GCAP + ! Save re-evaporation rate (kg/kg/s) + ! Convert from ( kg H2O/kg plume * mb plume ) -> ( kg H2O / kg air / s ) + kg_plume = 1d2 * ddraft * FMC1 * DXYP(j_debug) * bygrav + mc_dqevap(L,IC) = mc_dqevap(L,IC) + (dqevp/ddraft) * kg_plume * bykg(L) * bydtsrc +#endif + FSEVP = 0 if (PLK(L)*SMDN.gt.teeny) FSEVP = SLH*DQEVP/(PLK(L)*SMDN) SMDN=SMDN-SLH*DQEVP/PLK(L) @@ -2133,6 +2252,14 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) DTMOMR(:,L,1:NTX)=DTMOMR(:,L,1:NTX)-TMOM(:,L,1:NTX)*FENTRA #endif /* TRACERS_ON */ +#ifdef GCAP + !======================================================================== + ! Archive mass flux entrained to the downdraft (kg/m2/s) + !======================================================================== + mc_dd_ent( L, IC ) = 100.*EDRAFT*bygrav/dtsrc + mc_counts(L,IC,5) = mc_counts(L,IC,5) + 1 +#endif + else ! occasionally detrain into environment if ddraft too big FENTRA=EDRAFT/(DDRUP+teeny) ! < 0 DSM(L)=DSM(L)-FENTRA*SMDN @@ -2160,6 +2287,15 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) TMOMDN(xymoms,1:NTX)= TMOMDN(xymoms,1:NTX)*(1.+FENTRA) #endif /* TRACERS_ON */ +#ifdef GCAP + !======================================================================== + ! Archive mass flux detrained from the downdraft (kg/m2/s) + !======================================================================== + ! Note: Sign flipped + mc_dd_det( L, IC ) = -100.*EDRAFT*bygrav/dtsrc + mc_counts(L,IC,6) = mc_counts(L,IC,6) + 1 +#endif + end if end if @@ -2453,7 +2589,7 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) ! !**** reduce subsidence post hoc. ! LM1=max(1,L-1) ! if (TM(LM1,N)+TM(L,N).lt.0) then -! TM(L,N) = 0.d0 +! TM(L,N) = 0.d0 ! write(6,*) trname(n)," neg cannot be fixed!",L,TM(LM1:L,N) ! else ! TM(L-1,N)=TM(L-1,N)+TM(L,N) @@ -2469,7 +2605,7 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) enddo if(vsum.lt.0.) then #ifdef TRACERS_TOMAS - tm(l,n) = 0.d0 !dmw 2/1/2017 set neg tracer to zero + tm(l,n) = 0.d0 !dmw 2/1/2017 set neg tracer to zero write(6,*) trname(n)," neg cannot be fixed, setting to zero!",L,TM(1:L,N) #else write(6,*) trname(n)," neg cannot be fixed!",L,TM(1:L,N) @@ -2478,13 +2614,13 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) if(l-lborrow1.gt.1) then write(6,*) trname(n),' nonlocal borrow: it,i,j,l,tr,cm',itime,i_debug,j_debug,l,tm(lborrow1:l,n),cmneg(l) else - write(6,*) trname(n),' neg: it,i,j,l,tr,cm',itime,i_debug,j_debug,l,tm(l,n),cmneg(l) + write(6,*) trname(n),' neg: it,i,j,l,tr,cm',itime,i_debug,j_debug,l,tm(l,n),cmneg(l) #ifdef TRACERS_TOMAS - tm(l,n)=0.d0 !dmw 2/1/2017 set neg tracer to zero + tm(l,n)=0.d0 !dmw 2/1/2017 set neg tracer to zero ! setting tm=0 here means tm multiplier below is 1 write(6,*) 'TOMAS setting tm=0, not taking tm from below' -#endif - endif +#endif + endif ! note: borrowing from more than one layer is done by ! multiplication rather than subtraction tm(lborrow1:l-1,n)=tm(lborrow1:l-1,n)*(vsum/(vsum-tm(l,n))) @@ -2522,6 +2658,18 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) CUMFLX(L,IC,LMIN) = 100.*MCFLX(L)*bygrav/dtsrc DWNFLX(L,IC,LMIN) = 100.*DDMFLX(L)*bygrav/dtsrc endif +#endif +#ifdef GCAP + ! Archive plume mass fluxes (kg/m2/s) + + ! MCMFLX is the plume mass (CCM; mb) times fraction of plume area/gridbox area (FMC1) + mc_up_mf(L,IC) = mc_up_mf(L,IC) + 100.*MCFLX(L)*bygrav/dtsrc + mc_counts(L,IC,1) = mc_counts(L,IC,1) + 1 + + ! DDMFLX is the plume mass (DDM; mb) times fraction of plume area/gridbox area (FMC1) + mc_dd_mf(L,IC) = mc_dd_mf(L,IC) + 100.*DDMFLX(L)*bygrav/dtsrc + mc_counts(L,IC,4) = mc_counts(L,IC,4) + 1 + #endif end do @@ -2637,6 +2785,10 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) call stop_model("MSTCNV: negative cloud cover", 255) END IF +#ifdef GCAP + IF ( CCM(L) > 0 ) CCM_REF = CCM(L) +#endif + !**** PRECIPITATION IS ALLOWED TO EVAPORATE FULLY INTO A FRACTION OF !**** THE GRIDBOX HALF AS LARGE AS THE FRACTION OF GRIDBOX MASS THAT !**** CONVECTS @@ -2671,6 +2823,12 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) if (mcloud.gt.0) call get_dq_evap (smold(l),qmold(l),plk(l),airm(l),lhx,pl(l),prcp*AIRM(L)/MCLOUD, dqsum,fprcp) dqsum=dqsum*MCLOUD*BYAM(L) +#ifdef GCAP + ! Re-evaporated precipitation [kg/kg/s] + kg_plume = 1d2 * mcloud * fevap * DXYP(j_debug) * bygrav + if (mcloud.gt.0) mc_dqevap(L,IC) = mc_dqevap(L,IC) + (dqsum/mcloud) * kg_plume * bykg(L) * bydtsrc +#endif + PRCP=PRCP-DQSUM QM(L)=QM(L)+DQSUM @@ -2721,8 +2879,8 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) !Thus the grid-scale humidity is inversely weighted by the amount of !cloud in the grid box, and is assumed to be the average !of the humidity values before and after evaporation. - TNX1=(SM(L)*PLK(L)-SLH*DQSUM*(1./(2.*MCLOUD)-1.))*BYAM(L) - HEFF = Min (1d0, (QM(L)+DQSUM*(1/(2*MCLOUD)-1))*byAM(L) /QSAT(TNX1,LHX,PL(L))) + TNX1=(SM(L)*PLK(L)-SLH*DQSUM*(1./(2.*MCLOUD)-1.))*BYAM(L) + HEFF = Min (1d0, (QM(L)+DQSUM*(1/(2*MCLOUD)-1))*byAM(L) /QSAT(TNX1,LHX,PL(L))) Call GET_EVAP_FACTOR(NTX,TOLD,LHX,HEFF,FPRCP,FPRCPT,ntix) dtr(1:ntx) = fprcpt(1:ntx)*trprcp(1:ntx) @@ -2850,6 +3008,17 @@ subroutine MSTCNV(IERR,LERR,i_debug,j_debug) #endif #endif /* TRACERS_WATER */ +#ifdef GCAP + ! Convert (kg/kg)*mb to kg/m2/s + ! Use CCM_REF set to the last positive CCM(L) in the column to prevent + ! floating overflows for precipitation beneath the cloud base + IF ( TL(L) .gt. TF ) THEN + mc_pflx_l(L) = mc_pflx_l(L) + ( ( PRCP / CCM_REF ) * ( PL(L) * 1d2 / ( RGAS * TL(L) ) ) ) / dZm(L) + ELSE + mc_pflx_i(L) = mc_pflx_i(L) + ( ( PRCP / CCM_REF ) * ( PL(L) * 1d2 / ( RGAS * TL(L) ) ) ) / dZm(L) + ENDIF +#endif + !**** end of loop down from top of plume end do EVAP_PRECIP @@ -3104,7 +3273,7 @@ subroutine reset_tracer_work_arrays(l1,l2) end subroutine reset_tracer_work_arrays #endif /* TRACERS_ON */ - + end subroutine MSTCNV !**************************************************************************************** @@ -3328,6 +3497,7 @@ subroutine LSCOND(IERR,WMERR,LERR,i_debug,j_debug) integer K,L,N !@var K,L,N loop variables real*8 THBAR !@var THBAR potential temperature at layer edge + real*8 dqevap_tmp !**** !**** LARGE-SCALE CLOUDS AND PRECIPITATION @@ -3367,6 +3537,15 @@ subroutine LSCOND(IERR,WMERR,LERR,i_debug,j_debug) QHEATI=0. CLDSSL=0 TAUSSL=0 +#ifdef GCAP + CLDSSL3D=0. + TAUSSL3D=0. + ! Reset arrays that track the condensation/evaporation fluxes + ls_dqcond = 0d0 + ls_dqevap = 0d0 + ls_pflx_l(:) = 0d0 + ls_pflx_i(:) = 0d0 +#endif WMPR=0. prebar1=0. rh1=0. @@ -3877,11 +4056,15 @@ subroutine LSCOND(IERR,WMERR,LERR,i_debug,j_debug) !**** UNFAVORABLE CONDITIONS FOR CLOUDS TO EXIT, PRECIP OUT CLOUD WATER QHEATL(L)=0. ! QHEAT(L)=0. QHEATI(L)=0. + if (LHX.eq.LHE.and.QCLX(L).gt.0.) then call get_dq_evap (tl(l)*rh00(l)/plk(l),ql(l)*rh00(l),plk(l),rh00(l),lhx,pl(l),qclx(l)/(fssl(l)*rh00(l)), dqsum,fqcond1) DWDT=DQSUM*RH00(L)*FSSL(L) - +#ifdef GCAP + ls_dqevap(l) = ls_dqevap(l) + dwdt / dtsrc +#endif + !**** DWDT is amount of water going to vapour, store LH (sets QNEW below) QHEATL(L)=-DWDT*LHX*BYDTsrc if (debug) print*,"ls4",l,qheatl(l) @@ -3892,13 +4075,17 @@ subroutine LSCOND(IERR,WMERR,LERR,i_debug,j_debug) call get_dq_evap (tl(l)*rh00(l)/plk(l),ql(l)*rh00(l),plk(l),rh00(l),lhx,pl(l),qcix(l)/(fssl(l)*rh00(l)), dqsum,fqcond1) DWDT=DQSUM*RH00(L)*FSSL(L) - +#ifdef GCAP + ls_dqevap(l) = ls_dqevap(l) + dwdt / dtsrc +#endif + !**** DWDT is amount of water going to vapour, store LH (sets QNEW below) QHEATI(L)=-DWDT*LHX*BYDTsrc if (debug) print*,"ls4",l,qheat(l) PREP(L)=max(0d0,(QCIX(L)-DWDT)*BYDTsrc) ! precip out cloud water WMPR(L)=PREP(L)*DTsrc ! precip water (for opt. depth calculation) end if + ER(L)=(1.-RH(L))**ERP*LHX*PREBAR(L+1)*GbyAIRM0 ! GRAV/AIRM0 if(PREICE(L+1).gt.0..and.TL(L).lt.TF) & ER(L)=(1.-RHI)**ERP*LHX*PREBAR(L+1)*GbyAIRM0 ! GRAV/AIRM0 @@ -4298,7 +4485,10 @@ subroutine LSCOND(IERR,WMERR,LERR,i_debug,j_debug) SLH=LHX*BYSHA call get_dq_cond(tl(l),ql(l),1d0,1d0,lhx,pl(l),dqsum,fcond) - +#ifdef GCAP + ls_dqcond(l) = ls_dqcond(l) + dqsum * fssl(l) / dtsrc +#endif + if(DQSUM.gt.0.) then if (debug) print*,"lsB",l,slh*dqsum @@ -4501,9 +4691,17 @@ subroutine LSCOND(IERR,WMERR,LERR,i_debug,j_debug) HCNDSS=HCNDSS+FSSL(L)*(TL(L)-TOLD)*AIRM(L) SSHR(L)=SSHR(L)+FSSL(L)*(TL(L)-TOLD)*AIRM(L) DQLSC(L)=DQLSC(L)+FSSL(L)*(QL(L)-QOLD) + +#ifdef GCAP + ! Convert to kg/m2/s and from top edge (prebar) to bottom edge (ls_pflx*) + ls_pflx_l(l+1) = (prebar(l)-preice(l)) * 100. + ls_pflx_i(l+1) = preice(l) * 100. +#endif + end do CLOUD_FORMATION ! end of loop over L PRCPSS=max(0d0,PREBAR(1)*GRAV*DTsrc) ! fix small round off err + #ifdef TRACERS_WATER do n=1,ntx TRPRSS(n)=TRPRBAR(n,1) @@ -4982,6 +5180,19 @@ subroutine LSCOND(IERR,WMERR,LERR,i_debug,j_debug) #endif end do OPTICAL_THICKNESS +#ifdef GCAP + ! Archive these 3-D cloud fractions before converting to areal 2-D cloud fraction + DO L=1,LMCLD + IF ( TAUSSL(L) .lt. 0 ) THEN + CLDSSL3D(L) = 0d0 + TAUSSL3D(L) = 0d0 + ELSE + CLDSSL3D(L) = CLDSSL(L) + TAUSSL3D(L) = TAUSSL(L) + ENDIF + ENDDO +#endif + !**** CALCULATE OPTICAL THICKNESS do L=1,LMCLD CLDSV1(L)=CLDSSL(L) @@ -7067,7 +7278,7 @@ subroutine ISCCP_CLOUD_TYPES(sunlit,pfull & end subroutine get_dq_cond(sm,qm,plk,mass,lhx,pl, & ! input - dqsum,fcond) ! output + dqsum,fcond) !@sum get_dq calculation of condensation from vapour use CONSTANT, only : bysha implicit none @@ -7081,7 +7292,7 @@ subroutine get_dq_cond(sm,qm,plk,mass,lhx,pl, & ! input !@var DQSUM amount of water change (+ve is condensation) !@var FCOND fractional amount of vapour that condenses real*8, intent(OUT) :: dqsum,fcond - real*8 TP,QST,QSAT,DQSATDT,DQ,QMT,SLH + real*8 TP,QST,QSAT,DQSATDT,DQ,QMT,SLH,MP integer N DQSUM=0. diff --git a/model/CLOUDS2_DRV.F90 b/model/CLOUDS2_DRV.F90 index 6493f07..ae0baee 100644 --- a/model/CLOUDS2_DRV.F90 +++ b/model/CLOUDS2_DRV.F90 @@ -98,8 +98,9 @@ subroutine CONDSE use AMP_AEROSOL, only: NACTV #endif #endif -#if (defined CALCULATE_LIGHTNING) || (defined TRACERS_SPECIAL_Shindell) - USE LIGHTNING, only : FLASH_DENS, CG_DENS, FLASH_PERTURB,L440mbM1 +#if (defined CALCULATE_LIGHTNING) || (defined TRACERS_SPECIAL_Shindell) + USE LIGHTNING, only : FLASH_DENS, FLASH_PERTURB, L440mbM1 + USE LIGHTNING, only : CTH_SAVE #ifdef AUTOTUNE_LIGHTNING USE LIGHTNING, only : TUNE_LT_LAND, TUNE_LT_SEA USE LIGHTNING, only : LAND_FR_LIS, SEA_FR_LIS @@ -194,6 +195,10 @@ subroutine CONDSE ,wmclwp,wmctwp,CDNC_TOMAS #endif +#ifdef GCAP + use CLOUDS_COM, only : cldss3d + use CLOUDS, only : taussl3d, cldssl3d +#endif #if (defined CLD_AER_CDNC) || (defined CLD_SUBDD) use CLOUDS, only : cteml,cd3dl,cl3dl,ci3dl #endif @@ -203,6 +208,11 @@ subroutine CONDSE ! plume diagnostics use CLOUDS, only : CUMFLX,DWNFLX,WCUALL,ENTALL,DETALL, & MPLUMEALL,PLUME_MAX,PLUME_MIN +#endif +#ifdef GCAP + use CLOUDS, only : mc_up_mf, mc_dd_mf, mc_up_ent, mc_dd_ent, mc_up_det, mc_dd_det + use CLOUDS, only : mc_dqcond, mc_dqevap, ls_dqcond, ls_dqevap + use CLOUDS, only : mc_pflx_l, mc_pflx_i, ls_pflx_l, ls_pflx_i #endif use PBLCOM, only : dclev,egcm,w2gcm,pblht,pblptop use ATM_COM, only : pk,pek,pmid,pedn,gz,PMIDOLD,pdsig,MWs, & @@ -293,10 +303,30 @@ subroutine CONDSE integer LMIN #endif +#ifdef GCAP + real*8, & + dimension(GRID%I_STRT_HALO:GRID%I_STOP_HALO, & + GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM) :: & + dtrain, & ! Updraft detrainment flux in layer [kg/m2/s] + dqrcu, & ! Convective rainwater source in layer [kg/kg/s] + dqrlsan, & ! Stratiform rainwater source in layer [kg/kg/s] + reevapcn, & ! Evap./subl. of convective precip in layer [kg/kg/s] + reevapls ! Evap./subl. of stratiform precip in layer [kg/kg/s] + real*8, & + dimension(GRID%I_STRT_HALO:GRID%I_STOP_HALO, & + GRID%J_STRT_HALO:GRID%J_STOP_HALO,LM+1) :: & + cmfmc, & ! Upward cloud mass flux [kg/m2/s] + pficu, & ! Downward flux of convective ice precipitation [kg/m2/s] + pflcu, & ! Downward flux of convective liq precipitation [kg/m2/s] + pfilsan, & ! Downward flux of large-scale ice precipitation [kg/m2/s] + pfllsan ! Downward flux of large-scale liq precipitation [kg/m2/s] + integer LMIN +#endif + #ifdef CACHED_SUBDD #ifdef TRACERS_WATER !tracer precipitation variable name - character(len=20) :: trpname + character(len=20) :: trpname #endif #ifdef SCM ! isccp diagnostics save frequency histogram for subdd diagnostics @@ -541,8 +571,21 @@ subroutine CONDSE ! isccp frequency diags save_fq_isccp=0.d0 #endif -#endif - +#ifdef GCAP + ! plume diagnostics + dtrain = 0.d0 + dqrcu = 0.d0 + dqrlsan = 0.d0 + reevapcn = 0.d0 + reevapls = 0.d0 + cmfmc = 0.d0 + pficu = 0.d0 + pflcu = 0.d0 + pfilsan = 0.d0 + pfllsan = 0.d0 +#endif +#endif + call recalc_agrid_uv ! may not be necessary - check later ! @@ -625,7 +668,7 @@ subroutine CONDSE TUNE_LT_LAND, TUNE_LT_SEA !CALL STOP_MODEL( 'LTM Testing',17) FLASH_DENS = 0d0 - CG_DENS = 0d0 + !CG_DENS = 0d0 FLASH_UNC = 0d0 #endif @@ -890,7 +933,7 @@ subroutine CONDSE #if (defined CALCULATE_LIGHTNING) || (defined TRACERS_SPECIAL_Shindell) FLASH_DENS(i,j) = 0.d0 ! default for subdaily diag - CG_DENS(i,j) = 0.d0 ! default for subdaily diag + CTH_SAVE(i,j) = 0.d0 ! default for subdaily diag #ifdef AUTOTUNE_LIGHTNING FLASH_UNC(i,j) = 0.d0 ! default #endif @@ -1009,6 +1052,9 @@ subroutine CONDSE call inc_ajl(i,j,l,jl_mcdflx,DDMFLX(L)) call inc_ajl(i,j,l,jl_csizmc,CSIZEL(L)*CLDMCL(L)*AIRM(L)) aijl(i,j,l,ijl_MCamFX) = aijl(i,j,l,ijl_MCamFX) + MCFLX(L) +#ifdef GCAP + CMFMC(i,j,l+1) = MCFLX(l) * 100 * bygrav * bydtsrc ! mb -> kg m-2 s-1 +#endif end do do IT=1,NTYPE call INC_AJ(I,J,IT,J_PRCPMC,PRCPMC*FTYPE(IT,I,J)) @@ -1177,6 +1223,26 @@ subroutine CONDSE mc_pl_min_p2(I,J,:) = PLUME_MIN(2,:) endif #endif +#ifdef GCAP + ! plume diagnostics (kg/m2/s) + + !--------------------- + ! Mimic MERRA2 + !--------------------- + + ! Upward moist convective flux across top edge [kg/m2/s] + cmfmc(I,J,:) = ( mc_up_mf(:,1) + mc_up_mf(:,2) ) ! - & + ! ( mc_dd_mf(:,1) + mc_dd_mf(:,2) ) ! minus downdrafts + + ! Total detrainment flux in updrafts [kg/m2/s] + dtrain(I,J,:) = mc_up_det(:,1) + mc_up_det(:,2) + ! & + mc_dd_det(:,1) + mc_dd_det(:,2) ! including downdraft entrainment + + ! Condensation/evaporation rate from convective precipitation [kg/kg/s] + dqrcu(I,J,:) = mc_dqcond(:,1) + mc_dqcond(:,2)! Includes all cloud bases and both types + reevapcn(I,J,:) = mc_dqevap(:,1) + mc_dqevap(:,2) + +#endif #endif #ifdef TRACERS_ON @@ -1324,7 +1390,15 @@ subroutine CONDSE Itime,I,J,LERR,' CONDSE:H2O<0',WMERR,' ->0' !**** Accumulate diagnostics of LSCOND - +#ifdef GCAP + !--------------------- + ! Mimic MERRA2 + !--------------------- + ! Condensation/evaporation rate from large-scale precipitation [kg/kg/s] + dqrlsan(I,J,:) = ls_dqcond(:) + reevapls(I,J,:) = ls_dqevap(:) +#endif + #ifdef CLD_AER_CDNC ! code transplanted from LSCOND SMLWP=WMSUM @@ -1395,6 +1469,14 @@ subroutine CONDSE end do #endif +#ifdef GCAP + pflcu(i,j,:) = mc_pflx_l(:) ! kg m-2 s-1 + pfllsan(i,j,:) = ls_pflx_l(:) ! kg m-2 s-1 + pficu(i,j,:) = mc_pflx_i(:) ! kg m-2 s-1 + pfilsan(i,j,:) = ls_pflx_i(:) ! kg m-2 s-1 +#endif + + !**** CALCULATE PRECIPITATION HEAT FLUX (FALLS AT 0 DEGREES CENTIGRADE) !**** NEED TO TAKE ACCOUNT OF LATENT HEAT THOUGH if (LHP(1).ne.LHS) then @@ -1878,6 +1960,9 @@ subroutine CONDSE TAUSS(:,I,J)=TAUSSL(:) CLDSS(:,I,J)=CLDSSL(:) +#ifdef GCAP + CLDSS3D(:,I,J)=CLDSSL3D(:) +#endif CLDSAV(:,I,J)=CLDSAVL(:) CLDSAV1(:,I,J)=CLDSV1(:) SVLHX(:,I,J)=SVLHXL(:) @@ -2436,6 +2521,46 @@ subroutine CONDSE #ifdef CFMIP3_SUBDD case ('mc_lwp') call inc_subdd(subdd,k,cfmip_mc_lwp) +#endif +#if (defined CALCULATE_LIGHTNING) || (defined TRACERS_SPECIAL_Shindell) + case( 'FLASH_DENS' ) + do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr(i,j) = FLASH_DENS(i,j) + enddo; enddo + call inc_subdd(subdd,k,sddarr) + case( 'CTH' ) + do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr(i,j) = CTH_SAVE(i,j) + enddo; enddo + call inc_subdd(subdd,k,sddarr) +#endif +#ifdef GCAP + case('PRECANV') + sddarr(:,:) = 0 + call inc_subdd(subdd,k,sddarr) + case ('PRECCON') ! Normal mcp + do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr(i,j) = max(0.,prec(i,j)-precss(i,j))*bydtsrc ! Convert kg/m2 -> kg/m2/s + enddo; enddo + call inc_subdd(subdd,k,sddarr) + case( 'PRECLSC' ) + do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr(i,j) = precss(i,j)*bydtsrc ! Convert kg/m2 -> kg/m2/s + enddo; enddo + call inc_subdd(subdd,k,sddarr) + case( 'PRECSNO' ) + sddarr(:,:) = 0 + do j=j_0,j_1; do i=i_0,imaxj(j) + if(eprec(i,j).ge.0.) then + sddarr(i,j) = 0. + else + sddarr(i,j) = prec(i,j) + endif + enddo; enddo + call inc_subdd(subdd,k,sddarr) + case( 'PRECTOT' ) ! Based on prec + sddarr = prec * bydtsrc ! Convert kg/m2 -> kg/m2/s + call inc_subdd(subdd,k,sddarr) #endif end select enddo @@ -2457,6 +2582,12 @@ subroutine CONDSE #ifdef CFMIP3_SUBDD case ('mcamfx') call inc_subdd(subdd,k,cfmip_mcamfx) +#endif +#ifdef GCAP + case ('QL') + call inc_subdd(subdd,k,qcl) + case ('QI') + call inc_subdd(subdd,k,qci) #endif end select enddo @@ -2553,7 +2684,7 @@ subroutine CONDSE if(trpname.eq.trim(subdd%name(k))) then call inc_subdd(subdd,k,trprec(n,:,:)) exit ntm_loop - end if + end if end do ntm_loop end do !subdd diagnostics/variables end do !subdd groups @@ -2603,11 +2734,50 @@ subroutine CONDSE endif #endif +#ifdef GCAP + ! Plume diagnostics in aijl file + call find_groups('aijlh',grpids,ngroups) + do igrp=1,ngroups + subdd => subdd_groups(grpids(igrp)) + do k=1,subdd%ndiags + select case (subdd%name(k)) + case ('DTRAIN') + call inc_subdd(subdd,k, dtrain ) + case ('DQRCU') + call inc_subdd(subdd,k, dqrcu ) + case ('DQRLSAN') + call inc_subdd(subdd,k, dqrlsan ) + case ('REEVAPCN') + call inc_subdd(subdd,k, reevapcn ) + case ('REEVAPLS') + call inc_subdd(subdd,k, reevapls ) + end select + enddo + enddo + call find_groups('aijleh',grpids,ngroups) + do igrp=1,ngroups + subdd => subdd_groups(grpids(igrp)) + do k=1,subdd%ndiags + select case (subdd%name(k)) + case ('CMFMC') + call inc_subdd(subdd,k, cmfmc ) + case ('PFLCU') + call inc_subdd(subdd,k, pflcu ) + case ('PFLLSAN') + call inc_subdd(subdd,k, pfllsan ) + case ('PFICU') + call inc_subdd(subdd,k, pficu ) + case ('PFILSAN') + call inc_subdd(subdd,k, pfilsan ) + end select + enddo + enddo #endif - +#endif + #if (defined CALCULATE_LIGHTNING) || (defined TRACERS_SPECIAL_Shindell) #ifdef AUTOTUNE_LIGHTNING - CALL GLOBALSUM( grid, FLASH_UNC*(FOCEAN), SEA_FR_UNC(1), ALL=.true. ) + CALL GLOBALSUM( grid, FLASH_UNC*(FOCEAN), SEA_FR_UNC(1), ALL=.true. ) CALL GLOBALSUM( grid, FLASH_UNC*(1d0-FOCEAN), LAND_FR_UNC(1), ALL=.true. ) CNT_FR(1) = 1 #endif diff --git a/model/CLOUDS_COM.F90 b/model/CLOUDS_COM.F90 index 4abc8ac..5bb98cc 100644 --- a/model/CLOUDS_COM.F90 +++ b/model/CLOUDS_COM.F90 @@ -59,6 +59,10 @@ module CLOUDS_COM real*8, allocatable, dimension(:,:,:) :: TAUMC !@var CLDSS super-saturated cloud cover area (percent) real*8, allocatable, dimension(:,:,:) :: CLDSS +#ifdef GCAP +!@var CLDSS3D super-saturated cloud cover volume 3-D fraction + real*8, allocatable, dimension(:,:,:) :: CLDSS3D +#endif !@var CLDMC moist convective cloud cover area (percent) real*8, allocatable, dimension(:,:,:) :: CLDMC !@var CSIZMC,CSIZSS mc,ss effective cloud droplet radius (microns) @@ -147,7 +151,7 @@ subroutine get_cld_overlap (lmax, cldssl, cldmcl, CldTot, CldSS, CldMC, RandSS) ! overlap schemes are already realized by having picked a single MC ! random number for the whole column in addition to the ones for SS ! (That block was added to avoid unneeded computations in long runs) -! Note 2: Reverse looping over L was only kept for bit-wise consistency +! Note 2: Reverse looping over L was only kept for bit-wise consistency ! with the previous version of the code. if(present(RandSS)) then @@ -155,9 +159,9 @@ subroutine get_cld_overlap (lmax, cldssl, cldmcl, CldTot, CldSS, CldMC, RandSS) do L=lmax,1,-1 ! better: 1,lmax (and replace L+1 by L-1 below) if( cldssl(L) > 0 ) then if(same_cloud) RandSS(L) = RandSS(L+1) ! use same random # - same_cloud = .true. + same_cloud = .true. else - same_cloud = .false. + same_cloud = .false. end if end do end if @@ -168,14 +172,14 @@ subroutine get_cld_overlap (lmax, cldssl, cldmcl, CldTot, CldSS, CldMC, RandSS) if( cldssl(L) > 0 ) then !! if(same_cloud) RandSS(L) = RandSS(L-1) ! use same random # clearss_part = min( clearss_part, 1-cldssl(L) ) ! total overlap - same_cloud = .true. + same_cloud = .true. else ! clear sky layer - same_cloud = .false. + same_cloud = .false. clearss = clearss * clearss_part ! random overlap clearss_part = 1. ! reset for next cloud end if end do - clearss = clearss * clearss_part + clearss = clearss * clearss_part if( present(CldSS) ) CldSS = 1.-clearss !! Treat convective clouds as a single cloud with max. overlap @@ -219,6 +223,9 @@ subroutine ALLOC_CLOUDS_COM(grid) ULS,VLS,UMC,VMC,TLS,QLS,TAUSSIP,CSIZSSIP, & QLss,QIss,QLmc,QImc, & TMC,QMC,DDM1,AIRX,LMC,DDMS,TDN1,QDN1,DDML +#ifdef GCAP + use CLOUDS_COM, only : CLDSS3D +#endif #if (defined mjo_subdd) || (defined etc_subdd) use CLOUDS_COM, only : CLWC3D,CIWC3D,TLH3D,SLH3D,DLH3D,LLH3D #endif @@ -293,6 +300,10 @@ subroutine ALLOC_CLOUDS_COM(grid) QLmc(LM,I_0H:I_1H,J_0H:J_1H), & QImc(LM,I_0H:I_1H,J_0H:J_1H), & STAT=IER) +#ifdef GCAP + allocate( CLDSS3D(LM,I_0H:I_1H,J_0H:J_1H) ) + CLDSS3D = 0d0 +#endif #ifdef mjo_subdd allocate( & TMCDRY(LM,I_0H:I_1H,J_0H:J_1H), & diff --git a/model/GHY_COM.f b/model/GHY_COM.f index 688add7..1e20a60 100644 --- a/model/GHY_COM.f +++ b/model/GHY_COM.f @@ -136,6 +136,14 @@ MODULE GHY_COM !@var soil_surf_moist near surf soil moisture (kg/m^3) for subdd real*8, ALLOCATABLE, dimension(:,:) :: soil_surf_moist + +#ifdef GCAP +!@var Surface roughness height for diagnostics + REAL*8, allocatable, dimension(:,:) :: z0m_save +!@var LAI for diagnostics + REAL*8, allocatable, dimension(:,:) :: lai_save +#endif + END MODULE GHY_COM SUBROUTINE ALLOC_GHY_COM(grid) @@ -260,6 +268,13 @@ SUBROUTINE ALLOC_GHY_COM(grid) ALLOCATE( soil_surf_moist(I_0H:I_1H,J_0H:J_1H) ) soil_surf_moist(:,:) = 0.d0 +#ifdef GCAP + ALLOCATE( z0m_save(I_0H:I_1H,J_0H:J_1H) ) + z0m_save = 0d0 + ALLOCATE( lai_save(I_0H:I_1H,J_0H:J_1H) ) + lai_save = 0d0 +#endif + END SUBROUTINE ALLOC_GHY_COM subroutine read_landsurf_ic diff --git a/model/GHY_DRV.f b/model/GHY_DRV.f index dc3b5da..777c111 100644 --- a/model/GHY_DRV.f +++ b/model/GHY_DRV.f @@ -822,6 +822,9 @@ subroutine earth (ns,moddsf,moddd) #endif #ifdef TRACERS_GASEXCH_land_CO2 use tracer_com, only : n_CO2n +#endif +#ifdef GCAP + use ghy_com, only : wfcs, z0m_save, lai_save #endif use TimeConstants_mod, only: SECONDS_PER_YEAR implicit none @@ -1464,6 +1467,32 @@ subroutine earth (ns,moddsf,moddd) sddarr2d(i,j) = atmlnd%runo(i,j)*fearth(i,j) enddo; enddo call inc_subdd(subdd,k,sddarr2d) +#ifdef GCAP + case ('GWETTOP') + do j=j_0,j_1; do i=i_0,i_1 + sddarr2d(i,j)=0.0 + if (fearth(i,j).gt.0) then + sddarr2d(i,j)=(wearth(i,j)+aiearth(i,j))/(wfcs(i,j)+1e-20) + else + sddarr2d(i,j)=1 ! Set to 1 over oceans to match MERRA-2 + end if + enddo; enddo + call inc_subdd(subdd,k,sddarr2d) + case ('GWETROOT') + do j=j_0,j_1; do i=i_0,i_1 + sddarr2d(i,j)=0.0 + if (fearth(i,j).gt.0) then + sddarr2d(i,j)=(wearth(i,j)+aiearth(i,j))/(wfcs(i,j)+1e-20) + else + sddarr2d(i,j)=0 ! Set to 0 over oceans to match MERRA-2 + end if + enddo; enddo + call inc_subdd(subdd,k,sddarr2d) + case ('LAI') + call inc_subdd(subdd,k,lai_save) + case ('Z0M') + call inc_subdd(subdd,k,z0m_save) +#endif end select enddo enddo @@ -1649,6 +1678,9 @@ subroutine ghy_diag(i,j,jr,kr,ns,moddsf & ,airrig,aeirrig,alandC use ent_com, only : excess_C use ghy_com, only : gdeep, gsaveL, fearth +#ifdef GCAP + use ghy_com, only : lai_save +#endif USE CLOUDS_COM, only : DDMS use pbl_drv, only : t_pbl_args @@ -1759,6 +1791,9 @@ subroutine ghy_diag(i,j,jr,kr,ns,moddsf aij(i,j,ij_rauto)=aij(i,j,ij_rauto)+arauto*ptype aij(i,j,ij_clab)=aij(i,j,ij_clab)+(aclab/nisurf)*ptype aij(i,j,ij_lai)=aij(i,j,ij_lai)+(alai/nisurf)*ptype +#ifdef GCAP + lai_save(i,j) = (alai/nisurf)*ptype +#endif aij(i,j,ij_soilresp)=aij(i,j,ij_soilresp)+asoilresp*ptype aij(i,j,ij_soilCpoolsum)=aij(i,j,ij_soilCpoolsum) & + (asoilCpoolsum/nisurf)*ptype diff --git a/model/PBL_DRV.f b/model/PBL_DRV.f index 9589dc0..9d25732 100644 --- a/model/PBL_DRV.f +++ b/model/PBL_DRV.f @@ -68,7 +68,10 @@ SUBROUTINE PBL(I,J,IHC,ITYPE,PTYPE,pbl_args,atm) #ifdef SCM USE SCM_COM, only : SCMopt,SCMin #endif - +#ifdef GCAP + USE GHY_COM, only : z0m_save +#endif + use SOCPBL, only : npbl=>n, zgs, advanc USE PBLCOM @@ -244,7 +247,11 @@ SUBROUTINE PBL(I,J,IHC,ITYPE,PTYPE,pbl_args,atm) IF (ITYPE.GT.2) THEN Z0M=ROUGHL(I,J) ! 30./(10.**ROUGHL(I,J)) +#ifdef GCAP + z0m_save(I,J) = Z0M +#endif ENDIF + if ( isnan(zs1) ) zs1=0 ztop=zgs+zs1 ! zs1 is calculated before pbl is called c IF (pbl_args%TKV.EQ.pbl_args%TGV) c & pbl_args%TGV = 1.0001d0*pbl_args%TGV @@ -346,6 +353,9 @@ SUBROUTINE PBL(I,J,IHC,ITYPE,PTYPE,pbl_args,atm) & ,tr,trnradius,trndens,trnmm #endif & ) +#ifdef GCAP + z0m_save(i,j) = z0m +#endif atm%uabl(:,i,j)=upbl(:) atm%vabl(:,i,j)=vpbl(:) diff --git a/model/RADIATION.f b/model/RADIATION.f index 95e173d..15c6bf2 100644 --- a/model/RADIATION.f +++ b/model/RADIATION.f @@ -1774,6 +1774,9 @@ END SUBROUTINE RCOMPT SUBROUTINE RCOMPX use SURF_ALBEDO, only : getsur use O3mod, only : plbo3,nlo3,plbo3_traditional,nlo3_traditional +#ifdef GCAP + use O3mod, only : save_to3 +#endif #ifdef SCM use SCM_COM, only : SCMopt,SCMin #endif @@ -1858,6 +1861,10 @@ SUBROUTINE RCOMPX endif C-------------------------------- +#ifdef GCAP + ! Save DU of ozone used in calculation + save_to3(igcm,jgcm) = SUM(u0gas(:,3))*1000.0 +#endif C-------------------------------- if(set_aerosols_internally) then diff --git a/model/RAD_COM.f b/model/RAD_COM.f index 495d2f3..e018d76 100644 --- a/model/RAD_COM.f +++ b/model/RAD_COM.f @@ -93,6 +93,12 @@ MODULE RAD_COM !@var srnflb_save Net solar radiation (W/m^2) !@var trnflb_save Net thermal radiation (W/m^2) REAL*8,ALLOCATABLE,DIMENSION(:,:,:) :: srnflb_save,trnflb_save +#ifdef GCAP +!@var save_alb Surface albedo (unitless) + REAL*8,ALLOCATABLE,DIMENSION(:,:) :: save_alb +!@var TAUW3D,TAUI3D water,ice cloud opt. depths (for diags) + REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: TAUW3D,TAUI3D +#endif !@var TAUSUMW,TAUSUMI column-sum water,ice cloud opt. depths (for diags) REAL*8, DIMENSION(:,:), ALLOCATABLE :: TAUSUMW,TAUSUMI #ifdef mjo_subdd @@ -218,7 +224,7 @@ MODULE RAD_COM #endif !@var rad_to_chem save 3D quantities from radiation code for use in !@+ chemistry (or rest of model). 1=Ozone, 2=aerosol ext, 3=N2O, 4=CH4, -!@+ 5=CFC11+CFC12 +!@+ 5=CFC11+CFC12 REAL*8, ALLOCATABLE, DIMENSION(:,:,:,:) :: rad_to_chem !saved in rsf REAL*8, ALLOCATABLE, DIMENSION(:,:,:,:) :: rad_to_file #ifdef GCC_COUPLE_RAD @@ -392,6 +398,9 @@ SUBROUTINE ALLOC_RAD_COM(grid) #ifdef GCC_COUPLE_RAD * ,GCCco2_tracer_save,GCCco2rad_to_chem,GCCco2rad_to_file #endif +#ifdef GCAP + * ,save_alb,tauw3d,taui3d +#endif #ifdef mjo_subdd * ,SWHR_cnt,LWHR_cnt,SWHR,LWHR,OLR_acc,OLR_cnt * ,swu_avg,swu_cnt @@ -426,6 +435,10 @@ SUBROUTINE ALLOC_RAD_COM(grid) * DIFNIR(I_0H:I_1H, J_0H:J_1H), * TAUSUMW(I_0H:I_1H, J_0H:J_1H), * TAUSUMI(I_0H:I_1H, J_0H:J_1H), +#ifdef GCAP + * TAUW3D(I_0H:I_1H, J_0H:J_1H, LM), + * TAUI3D(I_0H:I_1H, J_0H:J_1H, LM), +#endif * SRDN(I_0H:I_1H, J_0H:J_1H), * CFRAC(I_0H:I_1H, J_0H:J_1H), * RCLD(LM, I_0H:I_1H, J_0H:J_1H), @@ -463,6 +476,13 @@ SUBROUTINE ALLOC_RAD_COM(grid) #endif * STAT=IER) +#ifdef GCAP +! Allocate and initialize array for holding surface albedo, which is only +! updated during daytime. + ALLOCATE( save_alb(I_0H:I_1H,J_0H:J_1H) ) + save_alb = 0. +#endif + #ifdef mjo_subdd OLR_acc=0. OLR_cnt=0. diff --git a/model/RAD_DRV.f b/model/RAD_DRV.f index 4e59e45..e843396 100644 --- a/model/RAD_DRV.f +++ b/model/RAD_DRV.f @@ -1641,6 +1641,9 @@ SUBROUTINE RADIA #ifdef GCC_COUPLE_RAD * ,GCCco2_tracer_save,GCCco2rad_to_chem #endif +#ifdef GCAP + * ,TAUW3D,TAUI3D +#endif #ifdef mjo_subdd * ,SWHR,LWHR,SWHR_cnt,LWHR_cnt,OLR_acc,OLR_cnt * ,swu_avg,swu_cnt @@ -1667,6 +1670,10 @@ SUBROUTINE RADIA * cldmc,cldss,csizmc,csizss,llow,lmid,lhi,fss,taussip,csizssip * ,QLss,QIss,QLmc,QImc * ,get_cld_overlap ! subroutine +#ifdef GCAP + use clouds_com, only : cldss3d + use constant, only : teeny +#endif USE DIAG_COM, only : ia_rad,jreg,aij=>aij_loc,aijl=>aijl_loc & ,ntype,ftype,itocean,itlake,itearth,itlandi,itoice,itlkice * ,adiurn=>adiurn_loc,ndiuvar,ia_rad_frc, @@ -1761,6 +1768,10 @@ SUBROUTINE RADIA #endif use DIAG_COM, only: ij_nintaerext,ij_nintaersca,ij_nintaerasy use RADPAR, only: nintaerext,nintaersca,nintaerasy +#ifdef GCAP + use RAD_COM, only : save_alb + use O3mod, only : save_to3 +#endif IMPLICIT NONE real*8 dz,rho C @@ -2239,6 +2250,11 @@ SUBROUTINE RADIA cfmip_qci = 0. cfmip_qcl = 0. #endif +#ifdef GCAP + tauw3d = 0. + taui3d = 0. +#endif + C**** C**** MAIN J LOOP C**** @@ -2399,6 +2415,9 @@ SUBROUTINE RADIA IF(SVLAT(L,I,J).EQ.LHE) THEN TAUWC(L)=cldx*TAUMCL OPTDW=OPTDW+TAUWC(L) +#ifdef GCAP + TAUW3D(I,J,L) = TAUW3D(I,J,L) + TAUMCL ! in-cloud vs. in-cell TAUWC(L) +#endif call inc_ajl(i,j,l,jl_wcld,1d0) call inc_ajl(i,j,l,jl_wcldwt,pdsig(l,i,j)) aij(i,j,ij_lwprad)=aij(i,j,ij_lwprad)+QLmc(l,i,j)*rhodz @@ -2415,6 +2434,9 @@ SUBROUTINE RADIA ELSE TAUIC(L)=cldx*TAUMCL OPTDI=OPTDI+TAUIC(L) +#ifdef GCAP + TAUI3D(I,J,L) = TAUI3D(I,J,L) + TAUMCL ! in-cloud vs. in-cell TAUIC(L) +#endif call inc_ajl(i,j,l,jl_icld,1d0) call inc_ajl(i,j,l,jl_icldwt,pdsig(l,i,j)) aij(i,j,ij_iwprad)=aij(i,j,ij_iwprad)+QImc(l,i,j)*rhodz @@ -2436,6 +2458,9 @@ SUBROUTINE RADIA IF(SVLHX(L,I,J).EQ.LHE) THEN TAUWC(L)=cldx*TAUSSL OPTDW=OPTDW+TAUWC(L) +#ifdef GCAP + TAUW3D(I,J,L) = TAUW3D(I,J,L) + TAUSSL ! in-cloud vs. in-cell TAUWC(L) +#endif call inc_ajl(i,j,l,jl_wcld,1d0) call inc_ajl(i,j,l,jl_wcldwt,pdsig(l,i,j)) aij(i,j,ij_lwprad)=aij(i,j,ij_lwprad)+QLss(l,i,j)*rhodz @@ -2453,6 +2478,9 @@ SUBROUTINE RADIA SIZEIC(L)=CSIZSSIP(L,I,J) TAUIC(L)=cldx*TAUSSLIP OPTDI=OPTDI+TAUIC(L) +#ifdef GCAP + TAUI3D(I,J,L) = TAUI3D(I,J,L) + TAUSSLIP ! in-cloud vs. in-cell TAUIC(L) +#endif call inc_ajl(i,j,l,jl_icld,1d0) call inc_ajl(i,j,l,jl_icldwt,pdsig(l,i,j)) aij(i,j,ij_iwprad)=aij(i,j,ij_iwprad)+QIss(l,i,j)*rhodz @@ -2469,6 +2497,9 @@ SUBROUTINE RADIA ELSE TAUIC(L)=cldx*TAUSSL OPTDI=OPTDI+TAUIC(L) +#ifdef GCAP + TAUI3D(I,J,L) = TAUI3D(I,J,L) + TAUSSL ! in-cloud vs. in-cell TAUIC(L) +#endif call inc_ajl(i,j,l,jl_icld,1d0) call inc_ajl(i,j,l,jl_icldwt,pdsig(l,i,j)) aij(i,j,ij_iwprad)=aij(i,j,ij_iwprad)+QIss(l,i,j)*rhodz @@ -2943,7 +2974,7 @@ SUBROUTINE RADIA end if #endif /* GCC_COUPLE_RAD */ #if (defined TRACERS_SPECIAL_Shindell) -! final (main) RCOMPX call can use tracer methane(or not): +! final (main) RCOMPX call can use tracer methane (or not): use_tracer_chem(2)=onoff_chem*Lmax_rad_CH4 if (is_set_param('initial_GHG_setup')) then call get_param('initial_GHG_setup', initial_GHG_setup) @@ -4258,6 +4289,74 @@ SUBROUTINE RADIA enddo enddo +#ifdef GCAP + call find_groups('aijlh',grpids,ngroups) + do igrp=1,ngroups + subdd => subdd_groups(grpids(igrp)) + do k=1,subdd%ndiags + select case (subdd%name(k)) + case ('OPTDEPTH') + do j=j_0,j_1; do i=i_0,imaxj(j); do l=1,lmaxsubdd + ! Weight mean cloud optical thickness by relative 2-D area fractions + sddarr3d(i,j,l)=(cldss(l,i,j)*TAUSS(l,i,j) + + & cldmc(l,i,j)*TAUMC(l,i,j) ) / + & ( cldss(l,i,j) + cldmc(l,i,j) + teeny ) + enddo; enddo; enddo + call inc_subdd(subdd,k,sddarr3d) + case ('CLOUD') + do j=j_0,j_1; do i=i_0,imaxj(j); do l=1,lmaxsubdd + sddarr3d(i,j,l)=min(1.0,CLDSS3D(l,i,j) + CLDMC(l,i,j)) + enddo; enddo; enddo + call inc_subdd(subdd,k,sddarr3d) + case ('TAUCLI') + call inc_subdd(subdd,k,taui3d) + case ('TAUCLW') + call inc_subdd(subdd,k,tauw3d) + end select + enddo + enddo + + call find_groups('aijh',grpids,ngroups) + do igrp=1,ngroups + subdd => subdd_groups(grpids(igrp)) + do k=1,subdd%ndiags + select case (subdd%name(k)) + case('PARDF') + do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr(i,j) = 0.82*srvissurf(i,j)*(1d0-fsrdir(i,j))*cosz1(i,j) + enddo; enddo + call inc_subdd(subdd,k,sddarr) + case('PARDR') + do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr(i,j) = 0.82*srvissurf(i,j)*(fsrdir(i,j))*cosz1(i,j) + enddo; enddo + call inc_subdd(subdd,k,sddarr) + case('ALBEDO') + do j=j_0,j_1; do i=i_0,imaxj(j) + if ( srdn(i,j).gt.0 ) then + ! Saves non-zero albedos for nighttime + save_alb(i,j) = 1d0-alb(i,j,1) + endif + sddarr(i,j) = save_alb(i,j) + enddo; enddo + call inc_subdd(subdd,k,sddarr) + case ('CLDTOT') ! totcld in standard model + call inc_subdd(subdd,k,cfrac) + case ('SWGDN') ! swds in standard model + do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr(i,j)=srdn(i,j)*cosz2(i,j) + enddo; enddo + call inc_subdd(subdd,k,sddarr) + case ('TO3') + do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr(i,j) = save_to3(i,j) + enddo; enddo + call inc_subdd(subdd,k,sddarr) + end select + enddo + enddo +#endif + call find_groups('rijlh',grpids,ngroups) do igrp=1,ngroups subdd => subdd_groups(grpids(igrp)) @@ -4266,15 +4365,15 @@ SUBROUTINE RADIA case ('MRCO2rad') do j=j_0,j_1; do i=i_0,imaxj(j); do l=1,lmaxsubdd sddarr3d(i,j,l) = CO2out(l,i,j) - enddo; enddo; enddo + enddo; enddo; enddo call inc_subdd(subdd,k,sddarr3d) case ('wtrtau') call inc_subdd(subdd,k,wtrtau) case ('icetau') call inc_subdd(subdd,k,icetau) - end select - end do - end do + end select + enddo + enddo #ifdef SCM call find_groups('rijlh',grpids,ngroups) @@ -5103,7 +5202,7 @@ subroutine rijh_defs(arr,nmax,decl_count) & units = 'W/m^2', & sched = sched_rad & ) -c +c arr(next()) = info_type_( & sname = 'swus', & lname = 'SOLAR UPWARD FLUX at SURFACE', diff --git a/model/RAD_UTILS.f b/model/RAD_UTILS.f index 297caca..3a6f204 100644 --- a/model/RAD_UTILS.f +++ b/model/RAD_UTILS.f @@ -3167,7 +3167,9 @@ module O3mod #ifdef HIGH_FREQUENCY_O3_INPUT type(timestream) :: OxHFstream,PSFforO3stream #endif - +#ifdef GCAP + REAL*8, ALLOCATABLE :: save_to3(:,:) +#endif !@dbparam use_sol_Ox_cycle if =1, a cycle of ozone is appled to !@+ o3year, as a function of the solar constant cycle. integer :: use_sol_Ox_cycle = 0 @@ -3263,6 +3265,12 @@ subroutine UPDO3D(JYEARO,JJDAYO,O3JDAY,O3JREF) plbo3(:)=plbo3_traditional(:) endif +#ifdef GCAP + allocate(save_to3(grid%i_strt:grid%i_stop, + & grid%j_strt:grid%j_stop)) + save_to3 = 0. +#endif + allocate(o3jday(nlo3,grid%i_strt:grid%i_stop, & grid%j_strt:grid%j_stop)) o3jday = 0. diff --git a/model/SEAICE_DRV.f b/model/SEAICE_DRV.f index fb41723..880dc8e 100644 --- a/model/SEAICE_DRV.f +++ b/model/SEAICE_DRV.f @@ -766,7 +766,7 @@ SUBROUTINE GROUND_SI(si_state,iceocn,atmice,atmocn) TRSNWIC = 0. #endif end if -c if (debug) write(6,'(A,2I4,4F11.6)') "si3",i,j,SNOW,1d3*SSIL(1) +c if (debug) write(6,'(A,2I4,4F11.6)') "si3",i,j,SNOW,1d3*SSIL(1) c $ /(XSI(1)*(ACE1I+SNOW)),1d3*SSIL(2)/(XSI(2)*(ACE1I+SNOW)) c $ ,1d3*(SSIL(1)+SSIL(2))/ACE1I c if (debug) print*,"si3",i,j,SNOW,HSIL,SSIL,MSI2, MSNWIC,HSNWIC, @@ -819,7 +819,7 @@ SUBROUTINE GROUND_SI(si_state,iceocn,atmice,atmocn) c * melt12,pond_melt(i,j),ti(HSIL(1)/(XSI(1)*(SNOW+ACE1I)), c * 1d3*SSIL(1)/(XSI(1)*(SNOW+ACE1I))), c * 1d3*(SSIL(1)+SSIL(2))/ACE1I - + C**** Net fluxes to ocean RUNOSI(I,J) = FMOC + RUN + MFLUX + MSNWIC ERUNOSI(I,J)= FHOC + ERUN + HFLUX + HSNWIC @@ -916,7 +916,7 @@ SUBROUTINE FORM_SI(si_state,iceocn,atmice) DO I=I_0,si_state%IMAXJ(J) c debug = i.eq.52.and.j.eq.87 - + PWATER=FWATER(I,J) ROICE=RSI(I,J) POICE=ROICE*PWATER @@ -1912,6 +1912,9 @@ SUBROUTINE READ_SEAICE(END_OF_DAY,atmocn,atmice) ZIMAX=2d0 IF(atmice%lat(i,j).GT.0.) ZIMAX=3.5d0 ! northern hemisphere rsinew = rsi(i,j) + if (rsi(i,j) > 1.d0) + * write (6,*) ' i,j, rsi(i,j) = ', i,j, rsi(i,j) + rsinew = min (rsi(i,j), 1.d0) MSINEW=RHOI*(ZIMIN-Z1I+(ZIMAX-ZIMIN)*RSINEW*DM(I,J)) C**** Ensure that lead fraction is consistent with kocean=1 case IF (RSINEW.gt.0) THEN @@ -1920,8 +1923,11 @@ SUBROUTINE READ_SEAICE(END_OF_DAY,atmocn,atmice) RSINEW = 1.-OPNOCN MSINEW=RHOI*(ZIMIN-Z1I+(ZIMAX-ZIMIN)*RSINEW*DM(I,J)) END IF - END IF + END IF RSI(I,J)=RSINEW + if (RSINEW > 1.d0) + * write (6,*) ' i,j, RSINEW = ', i,j, RSINEW + RSI(I,J)=min (RSINEW, 1.d0) MSI(I,J)=MSINEW END DO END DO diff --git a/model/SUBDD.f b/model/SUBDD.f index 3fc9b8b..84620bd 100644 --- a/model/SUBDD.f +++ b/model/SUBDD.f @@ -17,7 +17,7 @@ !@+ II. outputs whose registration is deferred until the stage !@+ of model execution at which the first time-slice of the !@+ output data is saved/accumulated. -!@+ +!@+ !@+ The type II interface is intended for outputs that do not !@+ naturally fall into one of the pre-existing categories of !@+ type I, and for which the creation of a new category @@ -31,7 +31,7 @@ !@+ parameters. Automatic vertical regridding of outputs to !@+ constant-pressure levels is currently only possible via I !@+ (will be added to II soon). -!@+ +!@+ !@+ Requests for type I outputs are made through rundeck !@+ strings SUBDD, SUBDD1, ..., following the traditional !@+ subdaily diagnostics framework. However, the parsing @@ -487,6 +487,8 @@ subroutine inc_subdd_3d(subdd,k,arr,jdim) endif call inc_subdd_alijh(subdd,k,arr) endif + case ('aijleh' ) + call inc_subdd_aijleh(subdd,k,arr) case ('aijph') if(jdim_ .eq. 2) then if(subdd%dsize3_input.eq.sizes(3)) then @@ -535,6 +537,30 @@ subroutine inc_subdd_aijlh(subdd,k,arr) return end subroutine inc_subdd_aijlh + subroutine inc_subdd_aijleh(subdd,k,arr) + use domain_decomp_atm, only : grid + use geom, only : imaxj + implicit none + type(subdd_type) :: subdd + integer :: k + real*8, dimension(grid%i_strt_halo:grid%i_stop_halo, + & grid%j_strt_halo:grid%j_stop_halo, + & subdd%dsize3_input) :: + & arr +c + integer :: i,j,l,ip +c + ip = subdd%subdd_period + do l=1,size(subdd%v5d,3) + do j=grid%j_strt,grid%j_stop + do i=grid%i_strt,imaxj(j) + subdd%v5d(i,j,l,ip,k) = subdd%v5d(i,j,l,ip,k) + arr(i,j,l) + enddo + enddo + enddo + return + end subroutine inc_subdd_aijleh + subroutine inc_subdd_alijh(subdd,k,arr) use domain_decomp_atm, only : grid use geom, only : imaxj @@ -1224,18 +1250,75 @@ subroutine create_group( & 'dist_im,dist_jm,lm_'//trim(grpname)// & ',nperiod_'//trim(grpname) endif +#ifdef GCAP +#ifdef CUBED_SPHERE + subdd%tile_dim_out = 4 + dimstr='(time,tile,lev,y,x) ;' +#else + dimstr='(time,lev,lat,lon) ;' +#endif +#else +#ifdef CUBED_SPHERE + subdd%tile_dim_out = 4 + dimstr='(time,tile,level,y,x) ;' +#else + dimstr='(time,level,lat,lon) ;' +#endif +#endif + allocate(lvlarr(dsize3)) + do l=1,dsize3 + lvlarr(l) = l + enddo +#ifdef GCAP + call add_coord(subdd%cdl0,'lev',size(lvlarr), + & long_name="levels", + & units = "1", + & coordvalues=lvlarr) +#else + call add_coord(subdd%cdl0,'level',size(lvlarr), + & coordvalues=lvlarr) +#endif + deallocate(lvlarr) + + case('aijleh' ) + if(catshape.eq.'aijleh') then + dsize3 = lmaxsubdd+1 + subdd%accshape = + & 'dist_im,dist_jm,lmaxsubdd+1,nperiod_'//trim(grpname) + else + dsize3 = dsize3_input + subdd%accshape = + & 'dist_im,dist_jm,lm_'//trim(grpname)// + & ',nperiod_'//trim(grpname) + endif +#ifdef GCAP +#ifdef CUBED_SPHERE + subdd%tile_dim_out = 4 + dimstr='(time,tile,lev,y,x) ;' +#else + dimstr='(time,lev,lat,lon) ;' +#endif +#else #ifdef CUBED_SPHERE subdd%tile_dim_out = 4 dimstr='(time,tile,level,y,x) ;' #else dimstr='(time,level,lat,lon) ;' +#endif #endif allocate(lvlarr(dsize3)) do l=1,dsize3 lvlarr(l) = l enddo +#ifdef GCAP + call add_coord(subdd%cdl0,'lev',size(lvlarr), + & long_name="levels", + & units = "1", + & coordvalues=lvlarr) +#else call add_coord(subdd%cdl0,'level',size(lvlarr), & coordvalues=lvlarr) +#endif deallocate(lvlarr) case('aijph') @@ -1367,7 +1450,7 @@ subroutine create_solo(subdd,nsubdd,vname,is_inst) allocate(subdd%name(ndiags)) allocate(subdd%denom(ndiags)) subdd%scale(:) = 1. - + subdd%name(1) = vname subdd%reduc(:) = reduc_avg subdd%denom(:) = 0 @@ -1376,7 +1459,7 @@ subroutine create_solo(subdd,nsubdd,vname,is_inst) subdd%is_inst = is_inst if(is_inst) then - subdd%sched(:) = sched_inst + subdd%sched(:) = sched_inst else subdd%sched(:) = sched_src endif @@ -1584,6 +1667,11 @@ subroutine parse_subdd input_sizes3(k) = lm call ijlh_defs(diaglists(1,k),nmax_possible,diaglens(k)) + k = k + 1 + catshapes(k) = 'aijleh'; categories(k) = 'aijleh' + input_sizes3(k) = lm+1 + call ijleh_defs(diaglists(1,k),nmax_possible,diaglens(k)) + k = k + 1 catshapes(k) = 'aijh'; categories(k) = 'rijh' input_sizes3(k) = 0 @@ -1846,7 +1934,7 @@ subroutine ijh_defs(arr,nmax,decl_count) & scale = real(nday,kind=8) & ) c - arr(next()) = info_type_( + arr(next()) = info_type_( & sname = 'evap', & lname = 'EVAPORATION', & units = 'mm/day', @@ -1982,7 +2070,7 @@ subroutine ijh_defs(arr,nmax,decl_count) & ) c arr(next()) = info_type_( - & sname = 'FLOPN', + & sname = 'FLOPN', & lname = 'Ice-Free Land Cover', & units = 'fraction' & ) @@ -2216,32 +2304,32 @@ subroutine ijh_defs(arr,nmax,decl_count) & sname = 'r_w_mc', & lname = 'Warm-Cloud effective Radius convective', & units = 'um' - & ) + & ) c arr(next()) = info_type_( & sname = 'r_i_mc', & lname = 'Ice-Cloud effective Radius convective', & units = 'um' - & ) + & ) c arr(next()) = info_type_( & sname = 'r_w_ls', & lname = 'Warm-Cloud effective Radius Large Scale', & units = 'um' - & ) + & ) c arr(next()) = info_type_( & sname = 'r_i_ls', & lname = 'Ice-Cloud effective Radius Large scale', & units = 'um' - & ) + & ) c arr(next()) = info_type_( & sname = 'pn', & lname = 'Number Concentration of dg > 0.1 um', & units = '#/m^2' & ) -c +c arr(next()) = info_type_( & sname = 'apn', & lname = 'Activated Particles Number Concentration', @@ -2259,6 +2347,255 @@ subroutine ijh_defs(arr,nmax,decl_count) & lname = 'Tropopause temperature', & units = 'K' & ) +c + arr(next()) = info_type_( + & sname = 'FLASH_DENS', + & lname = 'Lightning flash density', + & units = 'km-2 s-1' + & ) +c + arr(next()) = info_type_( + & sname = 'CTH', + & lname = 'Convective cloud depth', + & units = 'm' + & ) +c +#ifdef GCAP + arr(next()) = info_type_( + & sname = 'PS', + & lname = 'surface_pressure', + & units = 'Pa' + & ) + arr(next()) = info_type_( + & sname = 'FRSEAICE', + & lname = 'ice_covered_fraction_of_tile', + & units = '1' + & ) + arr(next()) = info_type_( + & sname = 'FRSNO', + & lname = 'fractional_area_of_land_snowcover', + & units = '1' + & ) + arr(next()) = info_type_( + & sname = 'GWETTOP', + & lname = 'surface_soil_wetness', + & units = '1' + & ) + arr(next()) = info_type_( + & sname = 'GWETROOT', + & lname = 'root_zone_soil_wetness', + & units = '1' + & ) + arr(next()) = info_type_( + & sname = 'HFLUX', + & lname = 'sensible_heat_flux_from_turbulence', + & units = 'W m-2' + & ) + arr(next()) = info_type_( + & sname = 'EFLUX', + & lname = 'total_latent_heat_flux', + & units = 'W m-2' + & ) + arr(next()) = info_type_( + & sname = 'LAI', + & lname = 'leaf_area_index', + & units = '1' + & ) + arr(next()) = info_type_( + & sname = 'LWI', + & lname = 'land_water_ice_flags', + & units = '1' + & ) + arr(next()) = info_type_( + & sname = 'PBLH', + & lname = 'planetary_boundary_layer_height', + & units = 'm' + & ) + arr(next()) = info_type_( + & sname = 'PRECANV', + & lname = 'anvil_precipitation', + & units = 'kg m-2 s-1' + & ) + arr(next()) = info_type_( + & sname = 'PRECCON', + & lname = 'convective_precipitation', + & units = 'kg m-2 s-1' + & ) + arr(next()) = info_type_( + & sname = 'PRECLSC', + & lname = 'nonanvil_large_scale_precipitation', + & units = 'kg m-2 s-1' + & ) + arr(next()) = info_type_( + & sname = 'PRECSNO', + & lname = 'snowfall', + & units = 'kg m-2 s-1' + & ) + arr(next()) = info_type_( + & sname = 'PRECTOT', + & lname = 'total_precipitation', + & units = 'kg m-2 s-1' + & ) + + arr(next()) = info_type_( + & sname = 'QV2M', + & lname = '2-meter_specific_humidity', + & units = 'kg kg-1' + & ) +! arr(next()) = info_type_( +! & sname = 'SEAICE00', +! & lname = 'sea_ice_area_fraction_0-10%', +! & units = '1' +! & ) +! arr(next()) = info_type_( +! & sname = 'SEAICE10', +! & lname = 'sea_ice_area_fraction_10-20%', +! & units = '1' +! & ) +! arr(next()) = info_type_( +! & sname = 'SEAICE20', +! & lname = 'sea_ice_area_fraction_20-30%', +! & units = '1' +! & ) +! arr(next()) = info_type_( +! & sname = 'SEAICE30', +! & lname = 'sea_ice_area_fraction_30-40%', +! & units = '1' +! & ) +! arr(next()) = info_type_( +! & sname = 'SEAICE40', +! & lname = 'sea_ice_area_fraction_40-50%', +! & units = '1' +! & ) +! arr(next()) = info_type_( +! & sname = 'SEAICE50', +! & lname = 'sea_ice_area_fraction_50-60%', +! & units = '1' +! & ) +! arr(next()) = info_type_( +! & sname = 'SEAICE60', +! & lname = 'sea_ice_area_fraction_60-70%', +! & units = '1' +! & ) +! arr(next()) = info_type_( +! & sname = 'SEAICE70', +! & lname = 'sea_ice_area_fraction_70-80%', +! & units = '1' +! & ) +! arr(next()) = info_type_( +! & sname = 'SEAICE80', +! & lname = 'sea_ice_area_fraction_80-90%', +! & units = '1' +! & ) +! arr(next()) = info_type_( +! & sname = 'SEAICE90', +! & lname = 'sea_ice_area_fraction_90-100%', +! & units = '1' +! & ) + arr(next()) = info_type_( + & sname = 'SLP', + & lname = 'sea_level_pressure', + & units = 'Pa' + & ) + arr(next()) = info_type_( + & sname = 'SNODP', + & lname = 'snow_depth', + & units = 'm' + & ) + arr(next()) = info_type_( + & sname = 'SNOMAS', + & lname = 'Total_snow_storage_land', + & units = 'kg m-2' + & ) + arr(next()) = info_type_( + & sname = 'TROPPT', + & lname = 'tropopause_pressure_based_on_thermal_estimate', + & units = 'Pa' + & ) + arr(next()) = info_type_( + & sname = 'TS', + & lname = 'surface_skin_temperature', + & units = 'K' + & ) + arr(next()) = info_type_( + & sname = 'T2M', + & lname = '2-meter_air_temperature', + & units = 'K' + & ) + arr(next()) = info_type_( + & sname = 'U10M', + & lname = '10-meter_eastward_wind', + & units = 'm s-1' + & ) + arr(next()) = info_type_( + & sname = 'USTAR', + & lname = 'surface_velocity_scale', + & units = 'm s-1' + & ) + arr(next()) = info_type_( + & sname = 'V10M', + & lname = '10-meter_northward_wind', + & units = 'm s-1' + & ) + arr(next()) = info_type_( + & sname = 'Z0M', + & lname = 'surface_roughness', + & units = 'm' + & ) + arr(next()) = info_type_( + & sname = 'T10M', + & lname = '10-meter_air_temperature', + & units = 'K' + & ) + arr(next()) = info_type_( + & sname = 'Q850', + & lname = 'specific_humidity_at_850_hPa', + & units = 'kg kg-1' + & ) + +! Rad scheduled diagnostics + arr(next()) = info_type_( + & sname = 'ALBEDO', + & lname = 'surface_albedo', + & units = '1', + & sched = sched_rad + & ) + arr(next()) = info_type_( + & sname = 'CLDTOT', + & lname = 'total_cloud_area_fraction', + & units = '1', + & sched = sched_rad + & ) +c +! EFLUX not used +c +! EVAP not used + arr(next()) = info_type_( + & sname = 'PARDF', + & lname = 'surface_downwelling_par_diffuse_flux', + & units = 'W/m^2', + & sched = sched_rad + & ) + arr(next()) = info_type_( + & sname = 'PARDR', + & lname = 'surface_downwelling_par_beam_flux', + & units = 'W/m^2', + & sched = sched_rad + & ) + arr(next()) = info_type_( + & sname = 'SWGDN', + & lname = 'surface_incoming_shortwave_flux', + & units = 'W m-2', + & sched = sched_rad + & ) + arr(next()) = info_type_( + & sname = 'TO3', + & lname = 'total_column_ozone', + & units = 'Dobsons', + & sched = sched_rad + & ) +#endif + return contains integer function next() @@ -2354,7 +2691,7 @@ subroutine ijph_defs(arr,nmax,decl_count) & lname = 'Raw model HDO', & units = 'kg/kg' & ) - + arr(next()) = info_type_( & sname = 'H2ORaw', & lname = 'Raw model H2O', @@ -2399,7 +2736,7 @@ subroutine ijph_defs(arr,nmax,decl_count) & dname = 'nTESGoodR', & units = 'kg/kg' & ) - + arr(next()) = info_type_( & sname = 'H2OR', & lname = 'H2O from retrieval-based TES operator', @@ -2649,12 +2986,103 @@ subroutine ijlh_defs(arr,nmax,decl_count) & units = 'cm-3' & ) #ifdef CFMIP3_SUBDD - arr(next()) = info_type_( + arr(next()) = info_type_( & sname = 'mcamfx', & lname = 'MC Air Mass Flux', & units = 'kg/s', & scale = 100.*bygrav/dtsrc & ) +#endif +#ifdef GCAP + arr(next()) = info_type_( + & sname = 'T', + & lname = 'air_temperature', + & units = 'K' + & ) + arr(next()) = info_type_( + & sname = 'QV', + & lname = 'specific_humidity', + & units = 'kg kg-1' + & ) + arr(next()) = info_type_( + & sname = 'DTRAIN', + & lname = 'detraining_mass_flux', + & units = 'kg m-s s-1' + & ) + arr(next()) = info_type_( + & sname = 'OMEGA', + & lname = 'vertical_pressure_velocity', + & units = 'Pa s-1' + & ) + arr(next()) = info_type_( + & sname = 'RH', + & lname = 'relative_humidity_after_moist', + & units = '1' + & ) + arr(next()) = info_type_( + & sname = 'U', + & lname = 'eastward_wind', + & units = 'm s-1' + & ) + arr(next()) = info_type_( + & sname = 'V', + & lname = 'northward_wind', + & units = 'm s-1' + & ) + arr(next()) = info_type_( + & sname = 'CLOUD', + & lname = 'cloud_fraction_for_radiation', + & units = '1', + & sched = sched_rad + & ) + arr(next()) = info_type_( + & sname = 'OPTDEPTH', + & lname = 'in_cloud_optical_thickness', + & units = '1', + & sched = sched_rad + & ) + arr(next()) = info_type_( + & sname = 'QI', + & lname = 'mass_fraction_of_cloud_ice_water', + & units = 'kg kg-1' + & ) + arr(next()) = info_type_( + & sname = 'QL', + & lname = 'mass_fraction_of_cloud_liquid_water', + & units = 'kg kg-1' + & ) + arr(next()) = info_type_( + & sname = 'TAUCLI', + & lname = 'in_cloud_optical_thickness_for_ice_clouds', + & units = '1', + & sched = sched_rad + & ) + arr(next()) = info_type_( + & sname = 'TAUCLW', + & lname = 'in_cloud_optical_thickness_for_liquid_clouds', + & units = '1', + & sched = sched_rad + & ) + arr(next()) = info_type_( + & sname = 'DQRCU', + & lname = 'convective_rainwater_source', + & units = 'kg kg-1 s-1' ! per grid-cell dry air + & ) + arr(next()) = info_type_( + & sname = 'DQRLSAN', + & lname = 'large_scale_rainwater_source', + & units = 'kg kg-1 s-1' ! per grid-cell dry air + & ) + arr(next()) = info_type_( + & sname = 'REEVAPCN', + & lname = 'evap_subl_of_convective_precipitation', + & units = 'kg kg-1 s-1' ! per grid-cell dry air + & ) + arr(next()) = info_type_( + & sname = 'REEVAPLS', + & lname = 'evap_subl_of_non_convective_precipitation', + & units = 'kg kg-1 s-1' ! per grid-cell dry air + & ) #endif return contains @@ -2664,6 +3092,58 @@ integer function next() end function next end subroutine ijlh_defs + subroutine ijleh_defs(arr,nmax,decl_count) +c +c 3D model-level edge outputs +c + use subdd_mod, only : info_type,sched_rad +! info_type_ is a homemade structure constructor for older compilers + use subdd_mod, only : info_type_ + use constant, only : bygrav,kapa + implicit none + integer :: nmax,decl_count + type(info_type) :: arr(nmax) +c +c note: next() is a locally declared function to increment decl_count +c + + decl_count = 0 + +#ifdef GCAP + arr(next()) = info_type_( + & sname = 'CMFMC', + & lname = 'cumulative_mass_flux', + & units = 'kg m-2 s-1' + & ) + arr(next()) = info_type_( + & sname = 'PFICU', + & lname = '3D_flux_of_ice_convective_precipitation', + & units = 'kg m-2 s-1' + & ) + arr(next()) = info_type_( + & sname = 'PFILSAN', + & lname = '3D_flux_of_ice_nonconvective_precipitation', + & units = 'kg m-2 s-1' + & ) + arr(next()) = info_type_( + & sname = 'PFLCU', + & lname = '3D_flux_of_liquid_convective_precipitation', + & units = 'kg m-2 s-1' + & ) + arr(next()) = info_type_( + & sname = 'PFLLSAN', + & lname = '3D_flux_of_liquid_nonconvective_precipitation', + & units = 'kg m-2 s-1' + & ) +#endif + + return + contains + integer function next() + decl_count = decl_count + 1 + next = decl_count + end function next + end subroutine ijleh_defs subroutine get_subdd_vinterp_coeffs use geom, only : imaxj @@ -3279,18 +3759,18 @@ subroutine read_subdd_rsf(fname) & rsf_equals_db_int('lmaxsubdd') & /) if(all(do_read)) then - do n=1,subdd_ngroups - grpname = subdd_groups(n)%grpname - call read_data(grid,fid,'nacc_'//trim(grpname), - & subdd_groups(n)%nacc, bcast_all=.true.) - if(allocated(subdd_groups(n)%v4d)) then - call read_dist_data(grid,fid,trim(grpname), - & subdd_groups(n)%v4d) - elseif(allocated(subdd_groups(n)%v5d)) then - call read_dist_data(grid,fid,trim(grpname), - & subdd_groups(n)%v5d) - endif - enddo + do n=1,subdd_ngroups + grpname = subdd_groups(n)%grpname + call read_data(grid,fid,'nacc_'//trim(grpname), + & subdd_groups(n)%nacc, bcast_all=.true.) + if(allocated(subdd_groups(n)%v4d)) then + call read_dist_data(grid,fid,trim(grpname), + & subdd_groups(n)%v4d) + elseif(allocated(subdd_groups(n)%v5d)) then + call read_dist_data(grid,fid,trim(grpname), + & subdd_groups(n)%v5d) + endif + enddo else if(am_i_root()) then write(6,*) 'WARNING: skipping rsf read of subdd info '// diff --git a/model/SURFACE.f b/model/SURFACE.f index 1dd4196..5449c33 100644 --- a/model/SURFACE.f +++ b/model/SURFACE.f @@ -1326,6 +1326,68 @@ SUBROUTINE SURFACE enddo; enddo call inc_subdd(subdd,k,sddarr2d) C +#ifdef GCAP + case ('FRSEAICE') + ! Based on FOICE + sddarr2d(:,:)=RSI(:,:)*FOCEAN(:,:) + call inc_subdd(subdd,k,sddarr2d) + case ('FRSNO') ! Fraction of snow cover, not including land ice sheets or ocean ice + do j=j_0,j_1; do i=i_0,imaxj(j) + xsubdd = 0. + ! Fraction of lake ice + if( snowi(i,j) > 0. ) + & xsubdd = xsubdd+rsi(i,j)*flake(i,j) + ! Fraction of land with snow cover + if( atmlnd%SNOWE(i,j) > 0. ) + & xsubdd = xsubdd + atmlnd%snowfr(i,j)*fearth(i,j) + xsubdd = min(1.0,xsubdd) + sddarr2d(i,j) = xsubdd + enddo; enddo + call inc_subdd(subdd,k,sddarr2d) + case ('LWI') + ! 0 = Ocean + Caspian Sea (in MERRA2, at least) + ! 1 = Land + Lakes + ! 2 = Ocean Ice + sddarr2d(:,:) = 1 + do j=j_0,j_1; do i=i_0,imaxj(j) + if ( focean(i,j) > fearth(i,j) ) sddarr2d(i,j) = 0 + if ( rsi(i,j)*focean(i,j) > 0.5 ) sddarr2d(i,j) = 2 + enddo; enddo + call inc_subdd(subdd,k,sddarr2d) + case ('TS') ! gtempr + do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr2d(i,j) = atmsrf%gtempr(i,j) + enddo; enddo + call inc_subdd(subdd,k,sddarr2d) + case ('T2M') ! tsavg + 273.15 + sddarr2d = atmsrf%tsavg(:,:)-tf+273.15 + call inc_subdd(subdd,k,sddarr2d) + case ('QV2M') + sddarr2d = atmsrf%qsavg(:,:) + call inc_subdd(subdd,k,sddarr2d) + case ('U10M') ! us + sddarr2d = atmsrf%usavg + call inc_subdd(subdd,k,sddarr2d) + case ('USTAR') ! ustar + do j=j_0,j_1; do i=i_0,imaxj(j) + sddarr2d(i,j) = atmsrf%ustar_pbl(i,j) + enddo; enddo + call inc_subdd(subdd,k,sddarr2d) + case ('V10M') ! vs + sddarr2d = atmsrf%vsavg + call inc_subdd(subdd,k,sddarr2d) + case ( 'SNODP' ) + sddarr2d(:,:) = 0. + ! atmsrf%SNOWDP = Snow thickness (m) over all layers + do j=j_0,j_1; do i=i_0,imaxj(j) + ! MERRA-2 has SNODP = 0 over ice sheets + sddarr2d(i,j) = atmsrf%SNOWDP(I,J) * ( 1d0 - flice(I,J) ) + enddo; enddo + call inc_subdd(subdd,k,sddarr2d) + case ( 'SNOMAS' ) ! Total snow storage over land (kg m-2) + sddarr2d = atmsrf%SNOW + call inc_subdd(subdd,k,sddarr2d) +#endif end select enddo enddo diff --git a/model/lightning.f b/model/lightning.f index 229d541..cbb70fe 100644 --- a/model/lightning.f +++ b/model/lightning.f @@ -32,6 +32,7 @@ module lightning ! SUBDD Instantaneous Diagnostic Arrays REAL*8, ALLOCATABLE, PUBLIC :: FLASH_DENS(:,:)! Flash density [fl/m2/s] REAL*8, ALLOCATABLE, PUBLIC :: CG_DENS(:,:) ! Cloud-to-ground density + REAL*8, ALLOCATABLE, PUBLIC :: CTH_SAVE(:,:) ! Cloud-top-height [km] #ifdef AUTOTUNE_LIGHTNING ! If AUTOTUNE_LIGHTNING is defined in the rundeck, then the tuning @@ -98,6 +99,7 @@ subroutine alloc_lightning(grid) USE ATM_COM, ONLY : PMIDL00 USE LIGHTNING, ONLY : CG_DENS USE LIGHTNING, ONLY : FLASH_DENS + USE LIGHTNING, ONLY : CTH_SAVE USE LIGHTNING, ONLY : L440mbM1 #ifdef TRACERS_SPECIAL_Shindell USE FILEMANAGER, ONLY : OPENUNIT, CLOSEUNIT @@ -131,6 +133,9 @@ subroutine alloc_lightning(grid) ! Allocate arrays for lightning !------------------------------------ + allocate( CTH_SAVE(I_0H:I_1H,J_0H:J_1H) ) + CTH_SAVE = 0d0 + allocate( FLASH_DENS(I_0H:I_1H,J_0H:J_1H) ) FLASH_DENS = 0d0 @@ -194,7 +199,7 @@ subroutine calc_lightning(i,j,lmax,lfrz,mflux,precon) USE LIGHTNING, ONLY : lightning_param USE LIGHTNING, ONLY : TUNE_LT_LAND, TUNE_LT_SEA, FLASH_PERTURB - USE LIGHTNING, ONLY : CG_DENS, FLASH_DENS, T_NEG_CTR + USE LIGHTNING, ONLY : CG_DENS, FLASH_DENS, T_NEG_CTR, CTH_SAVE #ifdef TRACERS_SPECIAL_Shindell USE LIGHTNING, ONLY : FLASH_YIELD_MIDLAT, FLASH_YIELD_TROPIC USE LIGHTNING, ONLY : CLDTOPL, ENOx_lgt @@ -468,8 +473,9 @@ subroutine calc_lightning(i,j,lmax,lfrz,mflux,precon) aij(i,j,ij_CtoG) = aij(i,j,ij_CtoG) + cg*DTsrc*byaxyp(i,j) ! Also save for SUBDD instantaneous output - FLASH_DENS(i,j) = flash*byaxyp(i,j) ! flashes/s -> flashes/m2/s - CG_DENS(i,j) = cg*byaxyp(i,j) ! flashes/s -> flashes/m2/s + FLASH_DENS(i,j) = flash*byaxyp(i,j)*1d6 ! flashes/s -> flashes/km2/s + CG_DENS(i,j) = cg*byaxyp(i,j)*1d6 ! flashes/s -> flashes/km2/s + CTH_SAVE(i,j) = htcon*1d3 ! km -> m #ifdef AUTOTUNE_LIGHTNING FLASH_UNC(i,j) = flashun ! flashes/s #endif diff --git a/templates/GCAP2.R b/templates/GCAP2.R new file mode 100644 index 0000000..0e9ceac --- /dev/null +++ b/templates/GCAP2.R @@ -0,0 +1,121 @@ +GCAP2.R GISS ModelE Lat-Lon Atmosphere Model, 1850 atm./ocean + +GCAP2 is based on E6F40 is based on LLF40 with updated aerosol/ozone input files for CMIP6 simulations + +Lat-lon: 2x2.5 degree horizontal resolution +F40: 40 vertical layers with standard hybrid coordinate, top at .1 mb +Atmospheric composition for year 1850 +Ocean climatology prescribed from years 1876-1885, CMIP6 +Uses turbulence scheme (no dry conv), grav.wave drag +Time steps: dynamics 3.75 min leap frog; physics 30 min.; radiation 2.5 hrs +Filters: U,V in E-W and N-S direction (after every physics time step) + U,V in E-W direction near poles (after every dynamics time step) + sea level pressure (after every physics time step) + +Preprocessor Options +#define STDHYB ! standard hybrid vertical coordinate +#define ATM_LAYERING L40 ! 40 layers, top at .1 mb +#define NEW_IO ! new I/O (netcdf) on +#define IRRIGATION_ON +#define MODIS_LAI +#define NEW_BCdalbsn +#define NEW_IO_SUBDD +#define CACHED_SUBDD +#define GCAP +#define CALCULATE_LIGHTNING +End Preprocessor Options + +Object modules: + ! resolution-specific source codes +Atm144x90 ! horizontal resolution is 144x90 -> 2x2.5deg +AtmLayering ! vertical resolution +DIAG_RES_F ! diagnostics +FFT144 ! Fast Fourier Transform + +IO_DRV ! new i/o + + ! GISS dynamics with gravity wave drag +ATMDYN MOMEN2ND ! atmospheric dynamics +QUS_DRV QUS3D ! advection of Q/tracers +STRATDYN STRAT_DIAG ! stratospheric dynamics (incl. gw drag) + +#include "latlon_source_files" +#include "modelE4_source_files" +#include "static_ocn_source_files" + +lightning +SUBDD + +Components: +#include "E4_components_nc" /* without "Ent" */ +Ent + +Component Options: +OPTS_Ent = ONLINE=YES PS_MODEL=FBB PFT_MODEL=ENT /* needed for "Ent" only */ +OPTS_dd2d = NC_IO=PNETCDF + +Data input files: +#include "IC_144x90_input_files" +#include "static_ocn_1880_144x90_input_files" +RVR=RD_Fd.nc ! river direction file +NAMERVR=RD_Fd.names.txt ! named river outlets + +#include "land144x90_input_files" +#include "rad_input_files" +#include "rad_144x90_input_files_CMIP6clim" + +MSU_wts=MSU_SSU_RSS_weights.txt ! MSU-diag +REG=REG2X2.5 ! special regions-diag + +Label and Namelist: (next 2 lines) +GCAP2 (LLF40 + updated aerosol/ozone input files for CMIP6 simulations, 1850 atm/ocean) + +&&PARAMETERS +#include "static_ocn_params" +#include "sdragF40_params" +#include "gwdragF40_params" + +! cond_scheme=2 ! newer conductance scheme (N. Kiang) ! not used with Ent + +! The following two lines are only used when aerosol/radiation interactions are off +FS8OPX=1.,1.,1.,1.,1.5,1.5,1.,1. +FT8OPX=1.,1.,1.,1.,1.,1.,1.3,1. + +! Increasing U00a decreases the high cloud cover; increasing U00b decreases net rad at TOA +U00a=0.655 ! above 850mb w/o MC region; tune this first to get 30-35% high clouds +U00b=1.00 ! below 850mb and MC regions; tune this last to get rad.balance +WMUI_multiplier = 2. +use_vmp=1 +radius_multiplier=1.1 + +PTLISO=0. ! pressure(mb) above which radiation assumes isothermal layers +H2ObyCH4=1. ! if =1. activates stratospheric H2O generated by CH4 without interactive chemistry +KSOLAR=2 ! 2: use long annual mean file ; 1: use short monthly file + +#include "atmCompos_1850_params" + +DTsrc=1800. ! cannot be changed after a run has been started +DT=225. +! parameters that control the Shapiro filter +DT_XUfilter=225. ! Shapiro filter on U in E-W direction; usually same as DT +DT_XVfilter=225. ! Shapiro filter on V in E-W direction; usually same as DT +DT_YVfilter=0. ! Shapiro filter on V in N-S direction +DT_YUfilter=0. ! Shapiro filter on U in N-S direction + +NIsurf=2 ! surface interaction computed NIsurf times per source time step +NRAD=1 ! radiation computed NRAD times per source time step +#include "diag_gcap2" + +Nssw=2 ! until diurnal diags are fixed, Nssw has to be even +Ndisk=960 ! write fort.1.nc or fort.2.nc every NDISK source time step +&&END_PARAMETERS + +&INPUTZ + YEARI=1949,MONTHI=12,DATEI=1,HOURI=0, ! pick IYEAR1=YEARI (default) or < YEARI + YEARE=1949,MONTHE=12,DATEE=2,HOURE=0, KDIAG=13*9, + ISTART=2,IRANDI=0, YEARE=1949,MONTHE=12,DATEE=1,HOURE=1, +/ +!! suggested settings for E6qsF40: +!! YEARI=1901,MONTHI=1,DATEI=1,HOURI=0, +!! YEARE=1931,MONTHE=1,DATEE=1,HOURE=0, KDIAG=12*0,9, +!! ISTART=8,IRANDI=0, YEARE=1901,MONTHE=1,DATEE=1,HOURE=1, diff --git a/templates/diag_gcap2 b/templates/diag_gcap2 new file mode 100644 index 0000000..437100a --- /dev/null +++ b/templates/diag_gcap2 @@ -0,0 +1,27 @@ +! parameters that affect at most diagn. output: standard if DTsrc=1800. (sec) +TAero_aod_diag=2 ! 0: no output; 1: save optical properties of TAero fields in aij for all bands; 2: save band6 only +aer_rad_forc=0 ! if set =1, radiation is called numerous times - slow !! +cloud_rad_forc=1 ! calls radiation twice; use =0 to save cpu time, 2= calculates crf_toa2 + +SUBDD='PS:6i QV:6i T:6i' +SUBDD1='ALBEDO:2 CLDTOT:2 EFLUX:2 FRSEAICE:2 FRSNO:2 GWETTOP:2' +SUBDD2='GWETROOT:2 HFLUX:2 LAI:2 LWI:2 PARDF:2 PARDR:2 PBLH:2' +SUBDD3='PRECANV:2 PRECCON:2 PRECLSC:2 PRECSNO:2 PRECTOT:2 SLP:2' +SUBDD4='SNODP:2 SNOMAS:2 SWGDN:2 T2M:2 TO3:2 TROPPT:2 TS:2' +SUBDD5='U10M:2 USTAR:2 V10M:2 Z0M:2 FLASH_DENS:2 CTH:2 QV2M:2' +SUBDD6='DTRAIN:6 OMEGA:6 RH:6 U:6 V:6' +SUBDD7='CLOUD:6 OPTDEPTH:6 QI:6 QL:6 TAUCLI:6 TAUCLW:6' +SUBDD8='DQRCU:6 DQRLSAN:6 REEVAPCN:6 REEVAPLS:6' +SUBDD9='CMFMC:6 PFICU:6 PFILSAN:6 PFLCU:6 PFLLSAN:6' + +NSUBDD=1 ! saving sub-daily diags every NSUBDD-th physics time step (1/2 hr) +DAYS_PER_FILE=1 + +KCOPY=1 ! 0: no output; 1: save .acc; 2: unused; 3: include ocean data +KRSF=12 ! 0: no output; X: save rsf at the beginning of every X month +isccp_diags=1 ! use =0 to save cpu time, but you lose some key diagnostics +nda5d=13 ! use =1 to get more accurate energy cons. diag (increases CPU time) +nda5s=13 ! use =1 to get more accurate energy cons. diag (increases CPU time) +ndaa=13 +nda5k=13 +nda4=48 ! to get daily energy history use nda4=24*3600/DTsrc