From 6443c5a204a39ffbb35f7c44f1ee93e85b1a3ce0 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Mon, 4 Dec 2023 12:04:42 -0500 Subject: [PATCH 1/3] Bypass MAPL_LoadBalance when not running load balancing --- GEOSsolar_GridComp/GEOS_SolarGridComp.F90 | 154 ++++++++++++---------- 1 file changed, 81 insertions(+), 73 deletions(-) diff --git a/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 b/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 index d3e6699..769c9f2 100644 --- a/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 +++ b/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 @@ -1433,7 +1433,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) real, pointer, dimension(:,:,:) :: ptr3d real, pointer, dimension(:,: ) :: ptr2d - + type (ESMF_State) :: AERO character(len=ESMF_MAXSTR) :: AS_FIELD_NAME integer :: AS_STATUS @@ -1497,10 +1497,10 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) logical :: PersistSolar logical :: do_no_aero_calc - + ! list of strings facility integer :: i - type S_ + type S_ character(len=:), allocatable :: str end type S_ type(S_), allocatable :: list(:) @@ -1610,7 +1610,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) ! Decide if should make OBIO exports !----------------------------------- - call MAPL_GetResource( MAPL, SOLAR_TO_OBIO, LABEL='SOLAR_TO_OBIO:', & + call MAPL_GetResource( MAPL, SOLAR_TO_OBIO, LABEL='SOLAR_TO_OBIO:', & DEFAULT=.FALSE., __RC__) ! Decide how to do solar forcing @@ -1740,10 +1740,10 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) ! Update the Sun position and weight the export variables ! ------------------------------------------------------- - if (UPDATE_FIRST) then - call MAPL_TimerOn (MAPL,"UPDATE",__RC__) - call UPDATE_EXPORT (IM,JM,LM,__RC__) - call MAPL_TimerOff (MAPL,"UPDATE",__RC__) + if (UPDATE_FIRST) then + call MAPL_TimerOn (MAPL,"UPDATE",__RC__) + call UPDATE_EXPORT (IM,JM,LM,__RC__) + call MAPL_TimerOff (MAPL,"UPDATE",__RC__) end if ! Periodically, refresh the internal state with a full solar calc @@ -1864,12 +1864,12 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) ! this line temporarily needed because of compiler bug allocate(list(1)); list(1) = S_('dummy') - + ! are without-aerosol exports requested? do_no_aero_calc = .false. list = [S_('FSWNA'), S_('FSWUNA'), S_('FSWDNA'), & S_('FSCNA'), S_('FSCUNA'), S_('FSCDNA'), & - S_('FSWBANDNA')] + S_('FSWBANDNA')] do i = 1, size(list) call MAPL_GetPointer( EXPORT, ptr3d, list(i)%str, __RC__) do_no_aero_calc = (do_no_aero_calc .or. associated(ptr3d)) @@ -2132,14 +2132,14 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) ! sub-gridscale condensate scaling for overlap scheme (ncols_block,nlay,ngpt) real(wp), dimension(:,:,:), allocatable :: zcw - + ! correlation length scales [m] for cloud presence and condensate (ncol) real, dimension(:), allocatable :: adl, rdl - + ! binomial probability of maximum overlap (cf. random overlap) ! for cloud presence and condensate (ncols_block,nlay-1) real(wp), dimension(:,:), allocatable :: alpha, rcorr - + ! TEMP ... see below real(wp) :: press_ref_min, ptop real(wp) :: temp_ref_min, tmin @@ -2236,7 +2236,10 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) !-------------------------------------------------------------------------- else - ZTH = max(.0001,ZTH) + ! This is the single-column model case. We *must* have lit points + if (adjustl(DYCORE)=="DATMO") then + ZTH = max(.0001,ZTH) + end if daytime = .true. NumLit = size(ZTH) end if @@ -2258,9 +2261,14 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) call MAPL_TimerOn(MAPL,"--CREATE") - call MAPL_BalanceCreate( & - OrgLen=NumLit, Comm=COMM, Handle=SolarBalanceHandle, & - BalLen=Num2do, BufLen=NumMax, __RC__) + if (LoadBalance) then + call MAPL_BalanceCreate( & + OrgLen=NumLit, Comm=COMM, Handle=SolarBalanceHandle, & + BalLen=Num2do, BufLen=NumMax, __RC__) + else + Num2do = NumLit + NumMax = NumLit + end if call MAPL_TimerOff(MAPL,"--CREATE") @@ -2270,15 +2278,15 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) ! column indicies Ig and Jg. ! The Outputs and InOuts are all INTERNAL variables. !-------------------------------------------------------------- - + NumImp = size(ImportSpec) NumInt = size(InternalSpec) - + ! Inputs to load balancing: ! All imports plus Ig, Jg, LATS, SLR & ZTH. ! Vertical only imports (PREF) are explicitly skipped later. NumInp = NumImp + 5 - + allocate( & SlicesInp(NumInp), NamesInp(NumInp), & SlicesInt(NumInt), NamesInt(NumInt), & @@ -2514,7 +2522,7 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) ! ----------------------- call MAPL_TimerOn(MAPL,"--DISTRIBUTE") - call MAPL_BalanceWork(BufInp,NumMax,Direction=MAPL_Distribute,Handle=SolarBalanceHandle,__RC__) + if (LoadBalance) call MAPL_BalanceWork(BufInp,NumMax,Direction=MAPL_Distribute,Handle=SolarBalanceHandle,__RC__) call MAPL_TimerOff(MAPL,"--DISTRIBUTE") ! @@@@@@@@@@@@@@@@@@@@@@ @@ -2525,17 +2533,17 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) ! the required length of their 1D buffers (BufInOut/BufOut). ! ---------------------------------------------------------- INT_VARS_1: do k=1,NumInt - - ! InOut or Out? + + ! InOut or Out? call MAPL_VarSpecGet(InternalSpec(k), & SHORT_NAME=short_name, DIMS=dims, UNGRIDDED_DIMS=ugdims, __RC__) ! later FAR variables will be InOut ... for now there are no InOut vars IntInOut(k) = .false. - ! save properties - NamesInt(k) = short_name + ! save properties + NamesInt(k) = short_name rgDim(k) = dims - + ! Skip vertical only variables. They dont require ! load-balancing since they have no horizontal dimension. if (dims == MAPL_DIMSVERTONLY) then @@ -2695,18 +2703,18 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) NIRR => ptr2(1:Num2do,1) case('DFNIRN') NIRF => ptr2(1:Num2do,1) - case('DRBANDN') + case('DRBANDN') DRBAND => ptr2(1:Num2do,:) - case('DFBANDN') + case('DFBANDN') DFBAND => ptr2(1:Num2do,:) - case('FSWNAN') - FSWA => ptr2(1:Num2do,:) - case('FSCNAN') - FSCA => ptr2(1:Num2do,:) - case('FSWUNAN') - FSWUA => ptr2(1:Num2do,:) - case('FSCUNAN') - FSCUA => ptr2(1:Num2do,:) + case('FSWNAN') + FSWA => ptr2(1:Num2do,:) + case('FSCNAN') + FSCA => ptr2(1:Num2do,:) + case('FSWUNAN') + FSWUA => ptr2(1:Num2do,:) + case('FSCUNAN') + FSCUA => ptr2(1:Num2do,:) case('FSWBANDNAN') FSWBANDA => ptr2(1:Num2do,:) case('COSZSW') @@ -2735,7 +2743,7 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) !---------------------------------- if (size(BufInOut) > 0) then call MAPL_TimerOn(MAPL,"--DISTRIBUTE") - call MAPL_BalanceWork(BufInOut,NumMax,Direction=MAPL_Distribute,Handle=SolarBalanceHandle,__RC__) + if (LoadBalance) call MAPL_BalanceWork(BufInOut,NumMax,Direction=MAPL_Distribute,Handle=SolarBalanceHandle,__RC__) call MAPL_TimerOff(MAPL,"--DISTRIBUTE") end if @@ -2746,7 +2754,7 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) call MAPL_TimerOn(MAPL,"-MISC") - ! report cosine solar zenith angle actually used by REFRESH + ! report cosine solar zenith angle actually used by REFRESH COSZSW = ZT ! save soon-to-be-calculated fluxes to correct set of internals @@ -2849,7 +2857,7 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) NIRR,NIRF,PARR,PARF,UVRR,UVRF, & FSWU,FSCU, & FSWBAND, & - SOLAR_TO_OBIO .and. include_aerosols, & + SOLAR_TO_OBIO .and. include_aerosols, & DRBAND, DFBAND, & __RC__ ) @@ -3101,10 +3109,10 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) ! cloud optics file is currently two-stream ! increment() will handle appropriate stream conversions allocate(ty_optical_props_2str::cloud_props_bnd,__STAT__) - + ! band-only initialization for pre-mcICA cloud optical properties TEST_(cloud_props_bnd%init(k_dist%get_band_lims_wavenumber())) - + ! g-point version for McICA sampled cloud optical properties select type (cloud_props_bnd) class is (ty_optical_props_2str) @@ -3118,12 +3126,12 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) call MAPL_GetResource( & MAPL, cloud_overlap_type, "RRTMGP_CLOUD_OVERLAP_TYPE_SW:", & DEFAULT='GEN_MAX_RAN_OVERLAP', __RC__) - + ! GEN_MAX_RAN_OVERLAP uses correlation lengths ! and possibly inhomogeneous condensate gen_mro = (cloud_overlap_type == "GEN_MAX_RAN_OVERLAP") if (gen_mro) then - + ! condensate inhomogeneous? ! see RadiationGC initialization cond_inhomo = condensate_inhomogeneous() @@ -3135,7 +3143,7 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) allocate(rdl(ncol),__STAT__) call correlation_length_condensate(ncol, ncol, doy, alat, rdl) endif - + endif ! =============================================================================== @@ -3200,14 +3208,14 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) ! set up aerosol optical properties need_aer_optical_props = (include_aerosols .and. implements_aerosol_optics) - if (need_aer_optical_props) then + if (need_aer_optical_props) then ! aerosol optics system is currently two-stream - ! increment() will handle appropriate stream conversions + ! increment() will handle appropriate stream conversions allocate(ty_optical_props_2str::aer_props,__STAT__) ! band-only initialization TEST_(aer_props%init(k_dist%get_band_lims_wavenumber())) end if - + !-------------------------------------------------------! ! Loop over blocks of blockSize columns ! ! - choose rrtmgp_blockSize for memory/time efficiency ! @@ -3349,12 +3357,12 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) if (need_aer_optical_props) then select type (aer_props) class is (ty_optical_props_2str) - + ! load un-normalized optical properties from aerosol system aer_props%tau = real(TAUA(colS:colE,:,:),kind=wp) aer_props%ssa = real(SSAA(colS:colE,:,:),kind=wp) aer_props%g = real(ASYA(colS:colE,:,:),kind=wp) - + ! renormalize where (aer_props%tau > 0._wp .and. aer_props%ssa > 0._wp) aer_props%g = aer_props%g / aer_props%ssa @@ -3364,7 +3372,7 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) aer_props%ssa = 0._wp aer_props%g = 0._wp end where - + class default TEST_('aerosol optical properties hardwired 2-stream for now') end select @@ -4032,13 +4040,13 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) if (size(BufOut) > 0) then call MAPL_TimerOn(MAPL,"--RETRIEVE") - call MAPL_BalanceWork(BufOut,NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) + if (LoadBalance) call MAPL_BalanceWork(BufOut,NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) call MAPL_TimerOff(MAPL,"--RETRIEVE") end if if (size(BufInOut) > 0) then call MAPL_TimerOn(MAPL,"--RETRIEVE") - call MAPL_BalanceWork(BufInOut,NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) + if (LoadBalance) call MAPL_BalanceWork(BufInOut,NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) call MAPL_TimerOff(MAPL,"--RETRIEVE") end if @@ -4047,14 +4055,14 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) ! resulting internals are then contiguous versions ! Note: InOut variables do not fill unmasked locations with a default, ! since the unmasked locations may contain potentially useful aged data. - + i1InOut = 1; i1Out = 1 INT_VARS_3: do k=1,NumInt if (SlicesInt(k) == 0) cycle - + if (IntInOut(k)) then pi1 => i1InOut - else + else pi1 => i1Out call MAPL_VarSpecGet(InternalSpec(k),DEFAULT=def,__RC__) endif @@ -4088,7 +4096,7 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) call ESMFL_StateGetPointerToData(INTERNAL,ptr3,NamesInt(k),__RC__) if (IntInOut(k)) then call UnPackIt(BufInOut(pi1),ptr3,daytime,NumMax,HorzDims,size(ptr3,3)) - else + else call UnPackIt(BufOut (pi1),ptr3,daytime,NumMax,HorzDims,size(ptr3,3),def) endif case(MAPL_DIMSHORZONLY) @@ -4110,7 +4118,7 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) deallocate(IntInOut,rgDim,ugDim,__STAT__) deallocate(BufInp,BufInOut,BufOut,__STAT__) call MAPL_TimerOn(MAPL,"--DESTROY") - call MAPL_BalanceDestroy(Handle=SolarBalanceHandle, __RC__) + if (LoadBalance) call MAPL_BalanceDestroy(Handle=SolarBalanceHandle, __RC__) call MAPL_TimerOff(MAPL,"--DESTROY") call MAPL_TimerOff(MAPL,"-BALANCE") @@ -5033,7 +5041,7 @@ subroutine UPDATE_EXPORT(IM,JM,LM, RC) if(associated(OSRCLR)) OSRCLR = (1.- FSCN(:,:, 0))*SLR if(associated( OSRNA)) OSRNA = (1.-FSWNAN(:,:, 0))*SLR if(associated(OSRCNA)) OSRCNA = (1.-FSCNAN(:,:, 0))*SLR - + ! SOLAR TO OBIO conversion ... ! Done in wavenum [cm^-1] space for reasons detailed under OBIO_bands_wavenum declaration ! --------------------------------------------------------------------------------------- @@ -5043,7 +5051,7 @@ subroutine UPDATE_EXPORT(IM,JM,LM, RC) ! first load SOLAR_bands_wavenum and specify ordering ! that makes it monotonically increasing ... - + if (USE_RRTMGP) then ! helper for testing RRTMGP error status on return; @@ -5084,7 +5092,7 @@ subroutine UPDATE_EXPORT(IM,JM,LM, RC) #undef TEST_ - elseif (USE_RRTMG) then + elseif (USE_RRTMG) then ! note: RRTMG bands are [wavenum1(jpb1:jpb2),wavenum2(jpb1:jpb2)] in [cm^-1] ! The index jpb1:jpb2 (16:29) is over the 14 bands. Band 14 is OUT of order. @@ -5210,19 +5218,19 @@ subroutine PackIt (Packed, UnPacked, MSK, Pdim, Udim, LM) logical, intent(IN ) :: MSK(Udim(1),Udim(2)) integer :: I, J, L, M - + do L = 1,LM - M = 1 - do J = 1,Udim(2) + M = 1 + do J = 1,Udim(2) do I = 1,Udim(1) if (MSK(I,J)) then Packed(M,L) = UnPacked(I,J,L) - M = M+1 - end if + M = M+1 + end if end do end do - end do - + end do + end subroutine PackIt ! Unpack masked locations from buffer @@ -5277,22 +5285,22 @@ subroutine choose_solar_scheme (MAPL, & USE_RRTMG = RFLAG /= 0. USE_CHOU = .not.USE_RRTMG end if - + _RETURN(_SUCCESS) - end subroutine choose_solar_scheme - + end subroutine choose_solar_scheme + subroutine choose_irrad_scheme (MAPL, & USE_RRTMGP, USE_RRTMG, USE_CHOU, & RC) - + type (MAPL_MetaComp), pointer, intent(in) :: MAPL logical, intent(out) :: USE_RRTMGP, USE_RRTMG, USE_CHOU integer, optional, intent(out) :: RC ! return code real :: RFLAG - integer :: STATUS - + integer :: STATUS + USE_RRTMGP = .false. USE_RRTMG = .false. USE_CHOU = .false. From d71d2645e562f1e8f94d2db68cfa1d8b27b9ce1d Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Mon, 4 Dec 2023 12:40:26 -0500 Subject: [PATCH 2/3] Add code to pass in maxpasses to the load balancer --- GEOSsolar_GridComp/GEOS_SolarGridComp.F90 | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 b/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 index 769c9f2..6c05b5b 100644 --- a/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 +++ b/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 @@ -1491,6 +1491,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) character(len=ESMF_MAXSTR) :: DYCORE integer :: SOLAR_LOAD_BALANCE integer :: SolarBalanceHandle + integer :: MaxPasses character(len=ESMF_MAXPATHLEN) :: SolCycFileName logical :: USE_NRLSSI2 @@ -1551,6 +1552,10 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) LoadBalance = .TRUE. end if + ! Note: We set the default to 100 as that is the default in MAPL_LoadBalance which + ! would have been used if not passed in + call MAPL_GetResource (MAPL, MaxPasses, 'SOLAR_LB_MAX_PASSES:', DEFAULT=100, __RC__) + ! Use time-varying co2 call ESMF_ClockGet(CLOCK, currTIME=CURRENTTIME, __RC__) call ESMF_TimeGet (CURRENTTIME, YY=YY, DayOfYear=DOY, __RC__) @@ -1891,6 +1896,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) call SORADCORE(IM,JM,LM, & include_aerosols = .false., & CURRTIME = CURRENTTIME+INTDT, & + MaxPasses = MaxPasses, & LoadBalance = LoadBalance, & __RC__) else @@ -1910,6 +1916,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) call SORADCORE(IM,JM,LM, & include_aerosols = .true., & CURRTIME = CURRENTTIME+INTDT,& + MaxPasses = MaxPasses, & LoadBalance = LoadBalance, & __RC__) @@ -1939,7 +1946,7 @@ subroutine RUN ( GC, IMPORT, EXPORT, CLOCK, RC ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) + subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,MaxPasses,LoadBalance,RC) ! RRTMGP module uses use mo_rte_kind, only: wp @@ -1979,6 +1986,7 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) integer, intent(IN ) :: IM, JM, LM logical, intent(IN ) :: include_aerosols type (ESMF_Time), intent(IN ) :: CURRTIME + integer, intent(IN ) :: MaxPasses logical, intent(IN ) :: LoadBalance integer, optional, intent(OUT) :: RC @@ -2263,7 +2271,7 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,LoadBalance,RC) if (LoadBalance) then call MAPL_BalanceCreate( & - OrgLen=NumLit, Comm=COMM, Handle=SolarBalanceHandle, & + OrgLen=NumLit, Comm=COMM, MaxPasses=MaxPasses, Handle=SolarBalanceHandle, & BalLen=Num2do, BufLen=NumMax, __RC__) else Num2do = NumLit From ef2ee2a6f0d23447f576233b6ea64f3003bca340 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Wed, 6 Dec 2023 12:39:08 -0500 Subject: [PATCH 3/3] Fix up MAPL profiler timers --- GEOSsolar_GridComp/GEOS_SolarGridComp.F90 | 69 +- .../rrtmg_sw/gcm_model/src/rrtmg_sw_rad.F90 | 721 +++++++++--------- 2 files changed, 389 insertions(+), 401 deletions(-) diff --git a/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 b/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 index 6c05b5b..4663f1b 100644 --- a/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 +++ b/GEOSsolar_GridComp/GEOS_SolarGridComp.F90 @@ -1313,40 +1313,6 @@ subroutine SetServices ( GC, RC ) !EOS - ! Set the Profiling timers - - call MAPL_TimerAdd(GC, name="PRELIMS" , __RC__) - call MAPL_TimerAdd(GC, name="REFRESH" , __RC__) - call MAPL_TimerAdd(GC, name="-AEROSOLS" , __RC__) - call MAPL_TimerAdd(GC, name="-SORAD" , __RC__) - call MAPL_TimerAdd(GC, name="--SORAD_RUN" , __RC__) - call MAPL_TimerAdd(GC, name="-RRTMG" , __RC__) - call MAPL_TimerAdd(GC, name="--RRTMG_RUN" , __RC__) - call MAPL_TimerAdd(GC, name="---RRTMG_PART" , __RC__) - call MAPL_TimerAdd(GC, name="---RRTMG_CLDSGEN" , __RC__) - call MAPL_TimerAdd(GC, name="---RRTMG_CLDPRMC" , __RC__) - call MAPL_TimerAdd(GC, name="---RRTMG_SETCOEF" , __RC__) - call MAPL_TimerAdd(GC, name="---RRTMG_TAUMOL" , __RC__) - call MAPL_TimerAdd(GC, name="---RRTMG_REFTRA" , __RC__) - call MAPL_TimerAdd(GC, name="---RRTMG_VRTQDR" , __RC__) - call MAPL_TimerAdd(GC, name="--RRTMG_INIT" , __RC__) - call MAPL_TimerAdd(GC, name="--RRTMG_FLIP" , __RC__) - call MAPL_TimerAdd(GC, name="-RRTMGP" , __RC__) - call MAPL_TimerAdd(GC, name="--RRTMGP_IO_GAS" , __RC__) - call MAPL_TimerAdd(GC, name="--RRTMGP_IO_CLOUDS" , __RC__) - call MAPL_TimerAdd(GC, name="--RRTMGP_CLOUD_OPTICS" , __RC__) - call MAPL_TimerAdd(GC, name="--RRTMGP_MCICA" , __RC__) - call MAPL_TimerAdd(GC, name="--RRTMGP_GAS_OPTICS" , __RC__) - call MAPL_TimerAdd(GC, name="--RRTMGP_RT" , __RC__) - call MAPL_TimerAdd(GC, name="--RRTMGP_POST" , __RC__) - call MAPL_TimerAdd(GC, name="-BALANCE" , __RC__) - call MAPL_TimerAdd(GC, name="--CREATE" , __RC__) - call MAPL_TimerAdd(GC, name="--DISTRIBUTE" , __RC__) - call MAPL_TimerAdd(GC, name="--RETRIEVE" , __RC__) - call MAPL_TimerAdd(GC, name="--DESTROY" , __RC__) - call MAPL_TimerAdd(GC, name="-MISC" , __RC__) - call MAPL_TimerAdd(GC, name="UPDATE" , __RC__) - ! Set Run method and use generic Initalize and Finalize methods call MAPL_GridCompSetEntryPoint (GC, ESMF_METHOD_RUN, Run, __RC__) call MAPL_GenericSetServices (GC, __RC__) @@ -2233,24 +2199,16 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,MaxPasses,LoadBalance,RC ! Identify lit soundings with the daytime mask !---------------------------------------------- - if (LoadBalance) then - daytime = ZTH > 0. - NumLit = count(daytime) - ! The load balancer does not work if there are no lit points. This is only ! important model-wise with the single-column model. Note we must protect ! ZTH since in solar, we divide by ZTH and, thus, we will get a divide-by- ! zero if not protected. !-------------------------------------------------------------------------- - else - ! This is the single-column model case. We *must* have lit points - if (adjustl(DYCORE)=="DATMO") then - ZTH = max(.0001,ZTH) - end if - daytime = .true. - NumLit = size(ZTH) - end if + if (adjustl(DYCORE)=="DATMO") ZTH = max(.0001,ZTH) + + daytime = ZTH > 0. + NumLit = count(daytime) ! Create a balancing strategy. This is a collective call on the communicator ! of the current VM. The original, unbalanced local work consists of (OrgLen) @@ -2749,11 +2707,11 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,MaxPasses,LoadBalance,RC ! Load balance the InOuts for Input !---------------------------------- + call MAPL_TimerOn(MAPL,"--DISTRIBUTE") if (size(BufInOut) > 0) then - call MAPL_TimerOn(MAPL,"--DISTRIBUTE") if (LoadBalance) call MAPL_BalanceWork(BufInOut,NumMax,Direction=MAPL_Distribute,Handle=SolarBalanceHandle,__RC__) - call MAPL_TimerOff(MAPL,"--DISTRIBUTE") end if + call MAPL_TimerOff(MAPL,"--DISTRIBUTE") call MAPL_TimerOff(MAPL,"-BALANCE") @@ -4046,17 +4004,12 @@ subroutine SORADCORE(IM,JM,LM,include_aerosols,CURRTIME,MaxPasses,LoadBalance,RC call MAPL_TimerOn(MAPL,"-BALANCE") - if (size(BufOut) > 0) then - call MAPL_TimerOn(MAPL,"--RETRIEVE") - if (LoadBalance) call MAPL_BalanceWork(BufOut,NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) - call MAPL_TimerOff(MAPL,"--RETRIEVE") - end if - - if (size(BufInOut) > 0) then - call MAPL_TimerOn(MAPL,"--RETRIEVE") - if (LoadBalance) call MAPL_BalanceWork(BufInOut,NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) - call MAPL_TimerOff(MAPL,"--RETRIEVE") + call MAPL_TimerOn(MAPL,"--RETRIEVE") + if (LoadBalance) then + if (size(BufOut) > 0) call MAPL_BalanceWork(BufOut, NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) + if (size(BufInOut) > 0) call MAPL_BalanceWork(BufInOut,NumMax,Direction=MAPL_Retrieve,Handle=SolarBalanceHandle,__RC__) end if + call MAPL_TimerOff(MAPL,"--RETRIEVE") ! Unpack the results. Fills "masked" (night) locations with default value from internal state !-------------------------------------------------------------------------------------------- diff --git a/GEOSsolar_GridComp/RRTMG/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.F90 b/GEOSsolar_GridComp/RRTMG/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.F90 index cf40cd9..e831f3a 100644 --- a/GEOSsolar_GridComp/RRTMG/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.F90 +++ b/GEOSsolar_GridComp/RRTMG/rrtmg_sw/gcm_model/src/rrtmg_sw_rad.F90 @@ -45,7 +45,7 @@ ! **************************************************************************** #include "MAPL_Generic.h" - + module rrtmg_sw_rad use ESMF @@ -92,7 +92,7 @@ subroutine rrtmg_sw (MAPL, & ! dimensions ! ---------- integer, intent(in) :: rpart ! Number of columns in a partition - integer, intent(in) :: ncol ! Number of horizontal columns + integer, intent(in) :: ncol ! Number of horizontal columns integer, intent(in) :: nlay ! Number of model layers ! orbit @@ -113,7 +113,7 @@ subroutine rrtmg_sw (MAPL, & ! ----------------- integer, intent(in) :: isolvar ! Flag for solar variability method - ! Notes: isolvar = -1 uses the Kurucz source function, while isolvar >= 0 use the NRLSSI2 + ! Notes: isolvar = -1 uses the Kurucz source function, while isolvar >= 0 use the NRLSSI2 ! solar model. First, the behavior for SCON = 0: isolvar in {-1,0,3} all have a fixed solar ! input at 1AU. For isolvar = -1 it is the Kurucz solar constant of 1368.22 Wm-2, while for ! isolvar = {0,3} it is the NRLSSI2 solar constant of 1360.85 Wm-2 (for the 100-50000 cm-1 @@ -152,7 +152,7 @@ subroutine rrtmg_sw (MAPL, & ! maximum (for 1=Mg, 2=SB respecively), and to vary linearly with solcycfrac between those ! extrema (see NRLSSI2 module for further details). ! Still discussing isolvar = 1, for SCON.eq.0 we take the hint from not explicitly - ! setting SCON to let an indsolvar.ne.1 choice cause a deviation from the internal solar + ! setting SCON to let an indsolvar.ne.1 choice cause a deviation from the internal solar ! constant, because, while the svar_{f,s} average to unity over a cycle *without* a time- ! varying indsolvar multiplier, they do not do so with it. If, on the other hand, a value ! SCON > 0 is provided, we ASSUME that it is a value that we should honor as a MEAN over @@ -272,7 +272,7 @@ subroutine rrtmg_sw (MAPL, & integer :: pncol integer :: STATUS ! for MAPL error reporting - + ! ASSERTs to catch unphysical or invalid inputs _ASSERT(all(play >= 0.), 'negative values in input: play') _ASSERT(all(plev >= 0.), 'negative values in input: plev') @@ -300,7 +300,7 @@ subroutine rrtmg_sw (MAPL, & else pncol = 2 end if - + ! do partitions call rrtmg_sw_sub (MAPL, & pncol, ncol, nlay, & @@ -319,9 +319,9 @@ subroutine rrtmg_sw (MAPL, & do_drfband, drband, dfband, & bndscl, indsolvar, solcycfrac, & ! optional inputs __RC__) - + _RETURN(_SUCCESS) - end subroutine rrtmg_sw + end subroutine rrtmg_sw subroutine rrtmg_sw_sub (MAPL, & @@ -391,24 +391,24 @@ subroutine rrtmg_sw_sub (MAPL, & ! cloud optics flags integer, intent(in) :: iceflgsw ! Flag for ice particle specifn integer, intent(in) :: liqflgsw ! Flag for liquid droplet specifn - + ! clouds real, intent(in) :: gcld (gncol,nlay) ! Cloud fraction real, intent(in) :: gciwp (gncol,nlay) ! In-cloud ice water path (g/m2) real, intent(in) :: gclwp (gncol,nlay) ! In-cloud liquid water path (g/m2) real, intent(in) :: grei (gncol,nlay) ! Cloud ice effective radius (um) real, intent(in) :: grel (gncol,nlay) ! Cloud drop effective radius (um) - + ! cloud overlap integer, intent(in) :: dyofyr ! Day of the year real, intent(in) :: gzm (gncol,nlay) ! Heights of level midpoints real, intent(in) :: galat (gncol) ! Latitudes of columns [radians] - + ! aerosols (optical props, non-delta-scaled) integer, intent(in) :: iaer ! aerosol flag (0=off, 10=on) - real, intent(in) :: gtauaer (gncol,nlay,nbndsw) ! aer optical depth (iaer=10 only) - real, intent(in) :: gssaaer (gncol,nlay,nbndsw) ! aer single scat alb (iaer=10 only) - real, intent(in) :: gasmaer (gncol,nlay,nbndsw) ! aer asymmetry param (iaer=10 only) + real, intent(in) :: gtauaer (gncol,nlay,nbndsw) ! aer optical depth (iaer=10 only) + real, intent(in) :: gssaaer (gncol,nlay,nbndsw) ! aer single scat alb (iaer=10 only) + real, intent(in) :: gasmaer (gncol,nlay,nbndsw) ! aer asymmetry param (iaer=10 only) ! surface albedos real, intent(in) :: gasdir (gncol) ! UV/vis surface albedo: direct rad @@ -416,7 +416,7 @@ subroutine rrtmg_sw_sub (MAPL, & real, intent(in) :: galdir (gncol) ! Near-IR surface albedo: direct rad real, intent(in) :: galdif (gncol) ! Near-IR surface albedo: diffuse rad - ! super-layer cloud fraction boundaries + ! super-layer cloud fraction boundaries integer, intent(in) :: cloudLM ! Low-mid integer, intent(in) :: cloudMH ! Mid-high @@ -473,12 +473,12 @@ subroutine rrtmg_sw_sub (MAPL, & ! surface albedos real :: albdir (nbndsw,pncol) ! surface albedo, direct real :: albdif (nbndsw,pncol) ! surface albedo, diffuse - + ! Atmosphere - setcoef ! -------------------- ! tropopause layer index - integer :: laytrop (pncol) + integer :: laytrop (pncol) ! gasesous absorbers real :: colh2o (nlay,pncol) ! column amount (h2o) @@ -489,17 +489,17 @@ subroutine rrtmg_sw_sub (MAPL, & real :: colmol (nlay,pncol) ! column amount (Rayleigh) ! continuum interpolation coefficients - integer :: indself (nlay,pncol) - integer :: indfor (nlay,pncol) - real :: selffac (nlay,pncol) - real :: selffrac (nlay,pncol) - real :: forfac (nlay,pncol) - real :: forfrac (nlay,pncol) + integer :: indself (nlay,pncol) + integer :: indfor (nlay,pncol) + real :: selffac (nlay,pncol) + real :: selffrac (nlay,pncol) + real :: forfac (nlay,pncol) + real :: forfrac (nlay,pncol) ! pressure and temperature interpolation coefficients integer, dimension (nlay,pncol) :: jp, jt, jt1 - real, dimension (nlay,pncol) :: fac00, fac01, fac10, fac11 - + real, dimension (nlay,pncol) :: fac00, fac01, fac10, fac11 + ! general real :: play (nlay, pncol) ! Layer pressures (hPa) real :: plev (nlay+1,pncol) ! Interface pressures (hPa) @@ -515,10 +515,10 @@ subroutine rrtmg_sw_sub (MAPL, & real :: clwp (nlay,pncol) ! In-cloud liq water path [g/m2] real :: rei (nlay,pncol) ! Cloud ice effective radius [um] real :: rel (nlay,pncol) ! Cloud drop effective radius [um] - + real :: alat (pncol) ! latitude for cloud overlap real :: zm (nlay,pncol) ! mid-layer hgt for cld overlap [m] - + logical :: cldymcl (nlay,ngptsw,pncol) ! cloud or not? [mcica] real :: ciwpmcl (nlay,ngptsw,pncol) ! in-cloud ice water path [mcica] [g/m2] real :: clwpmcl (nlay,ngptsw,pncol) ! in-cloud liq water path [mcica] [g/m2] @@ -528,7 +528,7 @@ subroutine rrtmg_sw_sub (MAPL, & real :: taormc (nlay,ngptsw,pncol) ! unscaled in-cloud optl depth [mcica] real :: ssacmc (nlay,ngptsw,pncol) ! in-cloud single scat albedo [mcica] real :: asmcmc (nlay,ngptsw,pncol) ! in-cloud asymmetry param [mcica] - + ! Atmosphere/clouds/aerosol - spcvrt,spcvmc ! ----------------------------------------- @@ -562,7 +562,7 @@ subroutine rrtmg_sw_sub (MAPL, & ! in-cloud PAR optical thicknesses real, dimension (pncol) :: ztautp, ztauhp, ztaump, ztaulp - + ! Solar variability multipliers ! ----------------------------- @@ -591,7 +591,7 @@ subroutine rrtmg_sw_sub (MAPL, & integer :: n, imol, gicol ! Loop indices real :: adjflx ! flux adjustment for Earth/Sun distance - + integer :: ipart, col_last, cols, cole, cc ! ncol is the actual number of gridcols in a partition, cf. pncol, @@ -601,7 +601,7 @@ subroutine rrtmg_sw_sub (MAPL, & ! other solar variability locals ! ------------------------------ real :: solvar (jpband) ! solar constant scaling factor by band - real :: indsolvar_scl (2) ! Adjusted facular and sunspot amplitude + real :: indsolvar_scl (2) ! Adjusted facular and sunspot amplitude ! scale factors (isolvar=1) real :: indsolvar_ndx (2) ! Facular and sunspot indices (isolvar=2) @@ -617,17 +617,17 @@ subroutine rrtmg_sw_sub (MAPL, & solvar(:) = 1. adjflux(:) = 1. svar_f = 1. - svar_s = 1. - svar_i = 1. - svar_f_bnd(:) = 1. - svar_s_bnd(:) = 1. - svar_i_bnd(:) = 1. + svar_s = 1. + svar_i = 1. + svar_f_bnd(:) = 1. + svar_s_bnd(:) = 1. + svar_i_bnd(:) = 1. ! isolvar == 1 specifies the position in AvgCyc11 through solcycfrac ! and allows scaling of solar cycle amplitudes as described in notes. ! ------------------------------------------------------------------ - if (isolvar .eq. 1) then + if (isolvar .eq. 1) then ! require solcycfrac present, else what's the point of using isolvar=1 ? if (.not.present(solcycfrac)) then @@ -635,10 +635,10 @@ subroutine rrtmg_sw_sub (MAPL, & end if solcycfr = solcycfrac - ! No amplitude scaling unless indsolvar is present. + ! No amplitude scaling unless indsolvar is present. indsolvar_scl(1:2) = 1. - if (present(indsolvar)) then + if (present(indsolvar)) then ! Adjust amplitude scaling of mean solar cycle to be unity at ! solar minimum (solcycfrac_min), to be the requested indsolvar @@ -654,15 +654,15 @@ subroutine rrtmg_sw_sub (MAPL, & ! isolvar == 2 allows direct specification of Mg and SB via indsolvar ! ------------------------------------------------------------------- - - if (isolvar .eq. 2) then + + if (isolvar .eq. 2) then ! default to mean indices indsolvar_ndx(1) = Mg_avg indsolvar_ndx(2) = SB_avg ! update to specified indices if provided - if (present(indsolvar)) then + if (present(indsolvar)) then indsolvar_ndx(1) = indsolvar(1) indsolvar_ndx(2) = indsolvar(2) endif @@ -675,7 +675,7 @@ subroutine rrtmg_sw_sub (MAPL, & ! Set flux adjustment for current Earth/Sun distance (two options) ! ---------------------------------------------------------------- - ! (Set adjflx to 1. to use constant Earth/Sun distance of 1 AU). + ! (Set adjflx to 1. to use constant Earth/Sun distance of 1 AU). ! 1) Provided by GCM via ADJES (from MAPL sun factor DIST ~ 1/r^2) adjflx = adjes @@ -692,7 +692,7 @@ subroutine rrtmg_sw_sub (MAPL, & ! and input solar constant SCON. ! -------------------------------------------------------- - if (scon == 0.) then + if (scon == 0.) then ! For scon = 0, use internally defined solar constant, which is ! 1368.22 Wm-2 (for ISOLVAR=-1) and 1360.85 Wm-2 (For ISOLVAR=0,3; @@ -718,7 +718,7 @@ subroutine rrtmg_sw_sub (MAPL, & elseif (isolvar .eq. 1) then ! Apply NRLSSI2 solar irradiance model at a specified solcycfr - ! within AvgCyc11, with the additional amplitude scalings in + ! within AvgCyc11, with the additional amplitude scalings in ! indsolvar_scl. ! interpolate mean solar cycle to solcycfr @@ -755,9 +755,9 @@ subroutine rrtmg_sw_sub (MAPL, & else _FAIL('invalid isolvar') - endif + endif - elseif (scon > 0.) then + elseif (scon > 0.) then ! Scale from internal to externally specified SCON. @@ -767,7 +767,7 @@ subroutine rrtmg_sw_sub (MAPL, & ! Scale from internal to requested solar constant. ! Apply optional scaling by band if bndscl present. - solvar(jpb1:jpb2) = scon / rrsw_scon + solvar(jpb1:jpb2) = scon / rrsw_scon if (present(bndscl)) & solvar(jpb1:jpb2) = solvar(jpb1:jpb2) * bndscl(:) @@ -775,7 +775,7 @@ subroutine rrtmg_sw_sub (MAPL, & ! Constant sun (NRLSSI2 model) ! Quiet sun, facular, and sunspot terms averaged over AvgCyc11. - ! Scale from internal to requested solar constant. + ! Scale from internal to requested solar constant. scon_int = Fint + Sint + Iint svar_r = scon / scon_int @@ -815,7 +815,7 @@ subroutine rrtmg_sw_sub (MAPL, & svar_f = (indsolvar_ndx(1) - Mg_0) / (Mg_avg - Mg_0) svar_s = (indsolvar_ndx(2) - SB_0) / (SB_avg - SB_0) - svar_i = (scon - (svar_f * Fint + svar_s * Sint)) / Iint + svar_i = (scon - (svar_f * Fint + svar_s * Sint)) / Iint elseif (isolvar .eq. 3) then @@ -835,7 +835,7 @@ subroutine rrtmg_sw_sub (MAPL, & else _FAIL('invalid isolvar') - endif + endif else _FAIL('scon cannot be negative!') @@ -849,7 +849,7 @@ subroutine rrtmg_sw_sub (MAPL, & if (isolvar < 0) then adjflux(jpb1:jpb2) = adjflux(jpb1:jpb2) * solvar(jpb1:jpb2) endif - + ! Build profile separation based on cloudiness, i.e., count and index ! clear/cloudy gridcolumns. The separation is based on whether the grid- ! column has cloud fraction in any layer (or not). This is based on the @@ -858,7 +858,7 @@ subroutine rrtmg_sw_sub (MAPL, & ! but the converse in not true ... can easily get a clear subcolumn for a ! cloudy gridcolumn. So, the gicol_clr can be assumed to yield all clear ! subcolumns, while the gicol_cld will yield both clear and cloudy sub- - ! columns. + ! columns. ncol_clr = 0 ncol_cld = 0 do gicol = 1,gncol @@ -887,7 +887,7 @@ subroutine rrtmg_sw_sub (MAPL, & do cc = 1,2 ! outer loop over clear then cloudy gridcolumns - if (cc == 1) then + if (cc == 1) then ! clear npart = npart_clr col_last = ncol_clr @@ -897,366 +897,401 @@ subroutine rrtmg_sw_sub (MAPL, & col_last = ncol_cld end if - ! loop over partitions - do ipart = 0,npart-1 + ! We need to call all the timeron/timeroff calls in the same order + ! so that all timers are called even if there are no lit points + ! (aka npart == 0) when the load balancer is off + if (npart == 0) then call MAPL_TimerOn(MAPL,"---RRTMG_PART",__RC__) + call MAPL_TimerOff(MAPL,"---RRTMG_PART",__RC__) - ! partition dimensions - cols = ipart * pncol + 1 - cole = (ipart + 1) * pncol - if (cole > col_last) cole = col_last - ncol = cole - cols + 1 + call MAPL_TimerOn(MAPL,"---RRTMG_CLDSGEN",__RC__) + call MAPL_TimerOff(MAPL,"---RRTMG_CLDSGEN",__RC__) - ! copy inputs into partition - ! -------------------------- + call MAPL_TimerOn(MAPL,"---RRTMG_CLDPRMC",__RC__) + call MAPL_TimerOff(MAPL,"---RRTMG_CLDPRMC",__RC__) - if (cc == 1) then + call MAPL_TimerOn(MAPL,"---RRTMG_SETCOEF",__RC__) + call MAPL_TimerOff(MAPL,"---RRTMG_SETCOEF",__RC__) - ! ----------------- - ! Clear gridcolumns - ! ----------------- + call MAPL_TimerOn(MAPL,"---RRTMG_TAUMOL",__RC__) + call MAPL_TimerOff(MAPL,"---RRTMG_TAUMOL",__RC__) - do icol = 1,ncol - gicol = gicol_clr(icol + cols - 1) + call MAPL_TimerOn(MAPL,"---RRTMG_REFTRA",__RC__) + call MAPL_TimerOff(MAPL,"---RRTMG_REFTRA",__RC__) + + call MAPL_TimerOn(MAPL,"---RRTMG_VRTQDR",__RC__) + call MAPL_TimerOff(MAPL,"---RRTMG_VRTQDR",__RC__) + else - ! assign surface albedos to bands + ! loop over partitions + do ipart = 0,npart-1 - ! near IR bands 14=nbndsw and 1-8 - ! 820-12850 cm-1, 0.778-12.2 um - do ibnd=1,8 - albdir(ibnd,icol) = galdir(gicol) - albdif(ibnd,icol) = galdif(gicol) - enddo - albdir(nbndsw,icol) = galdir(gicol) - albdif(nbndsw,icol) = galdif(gicol) - - ! UV/Vis bands 10-13 - ! 16000-50000 cm-1, 0.200-0.625 um - do ibnd=10,13 - albdir(ibnd,icol) = gasdir(gicol) - albdif(ibnd,icol) = gasdif(gicol) - enddo + call MAPL_TimerOn(MAPL,"---RRTMG_PART",__RC__) - ! Transition band 9 - ! 12850-16000 cm-1, 0.625-0.778 um - ! Take average, dmlee - albdir(9,icol) = (gasdir(gicol)+galdir(gicol))/2. - albdif(9,icol) = (gasdif(gicol)+galdif(gicol))/2. + ! partition dimensions + cols = ipart * pncol + 1 + cole = (ipart + 1) * pncol + if (cole > col_last) cole = col_last + ncol = cole - cols + 1 - enddo + ! copy inputs into partition + ! -------------------------- - ! copy in partition (general) - do icol = 1,ncol - gicol = gicol_clr(icol + cols - 1) - - play(:,icol) = gplay(gicol,1:nlay) - plev(:,icol) = gplev(gicol,1:nlay+1) - tlay(:,icol) = gtlay(gicol,1:nlay) - coszen(icol) = gcoszen(gicol) + if (cc == 1) then - enddo + ! ----------------- + ! Clear gridcolumns + ! ----------------- - ! copy in partition (aerosols) - if (iaer == 10) then do icol = 1,ncol gicol = gicol_clr(icol + cols - 1) - do ibnd = 1,nbndsw - taua(1:nlay,ibnd,icol) = gtauaer(gicol,1:nlay,ibnd) - asya(1:nlay,ibnd,icol) = gasmaer(gicol,1:nlay,ibnd) - omga(1:nlay,ibnd,icol) = gssaaer(gicol,1:nlay,ibnd) + + ! assign surface albedos to bands + + ! near IR bands 14=nbndsw and 1-8 + ! 820-12850 cm-1, 0.778-12.2 um + do ibnd=1,8 + albdir(ibnd,icol) = galdir(gicol) + albdif(ibnd,icol) = galdif(gicol) + enddo + albdir(nbndsw,icol) = galdir(gicol) + albdif(nbndsw,icol) = galdif(gicol) + + ! UV/Vis bands 10-13 + ! 16000-50000 cm-1, 0.200-0.625 um + do ibnd=10,13 + albdir(ibnd,icol) = gasdir(gicol) + albdif(ibnd,icol) = gasdif(gicol) enddo - enddo - endif - ! copy in partition (gases) - do icol = 1,ncol - gicol = gicol_clr(icol + cols - 1) - colh2o(:,icol) = gh2ovmr(gicol,1:nlay) - colco2(:,icol) = gco2vmr(gicol,1:nlay) - colo3 (:,icol) = go3vmr (gicol,1:nlay) - colch4(:,icol) = gch4vmr(gicol,1:nlay) - colo2 (:,icol) = go2vmr (gicol,1:nlay) - end do - - else - - ! ------------------ - ! Cloudy gridcolumns - ! ------------------ - - do icol = 1,ncol - gicol = gicol_cld(icol + cols - 1) - - ! assign surface albedos to bands - - ! near IR bands 14=nbndsw and 1-8 - ! 820-12850 cm-1, 0.778-12.2 um - do ibnd=1,8 - albdir(ibnd,icol) = galdir(gicol) - albdif(ibnd,icol) = galdif(gicol) + ! Transition band 9 + ! 12850-16000 cm-1, 0.625-0.778 um + ! Take average, dmlee + albdir(9,icol) = (gasdir(gicol)+galdir(gicol))/2. + albdif(9,icol) = (gasdif(gicol)+galdif(gicol))/2. + enddo - albdir(nbndsw,icol) = galdir(gicol) - albdif(nbndsw,icol) = galdif(gicol) - - ! UV/Vis bands 10-13 - ! 16000-50000 cm-1, 0.200-0.625 um - do ibnd=10,13 - albdir(ibnd,icol) = gasdir(gicol) - albdif(ibnd,icol) = gasdif(gicol) + + ! copy in partition (general) + do icol = 1,ncol + gicol = gicol_clr(icol + cols - 1) + + play(:,icol) = gplay(gicol,1:nlay) + plev(:,icol) = gplev(gicol,1:nlay+1) + tlay(:,icol) = gtlay(gicol,1:nlay) + coszen(icol) = gcoszen(gicol) + enddo - ! Transition band 9 - ! 12850-16000 cm-1, 0.625-0.778 um - ! Take average, dmlee - albdir(9,icol) = (gasdir(gicol)+galdir(gicol))/2. - albdif(9,icol) = (gasdif(gicol)+galdif(gicol))/2. + ! copy in partition (aerosols) + if (iaer == 10) then + do icol = 1,ncol + gicol = gicol_clr(icol + cols - 1) + do ibnd = 1,nbndsw + taua(1:nlay,ibnd,icol) = gtauaer(gicol,1:nlay,ibnd) + asya(1:nlay,ibnd,icol) = gasmaer(gicol,1:nlay,ibnd) + omga(1:nlay,ibnd,icol) = gssaaer(gicol,1:nlay,ibnd) + enddo + enddo + endif - enddo - - ! copy in partition (general and cloud physical props) - do icol = 1,ncol - gicol = gicol_cld(icol + cols - 1) - - play(:,icol) = gplay(gicol,1:nlay) - plev(:,icol) = gplev(gicol,1:nlay+1) - tlay(:,icol) = gtlay(gicol,1:nlay) - cld (:,icol) = gcld (gicol,1:nlay) - ciwp(:,icol) = gciwp(gicol,1:nlay) - clwp(:,icol) = gclwp(gicol,1:nlay) - rei (:,icol) = grei (gicol,1:nlay) - rel (:,icol) = grel (gicol,1:nlay) - zm (:,icol) = gzm (gicol,1:nlay) - alat (icol) = galat (gicol) - coszen(icol) = gcoszen(gicol) - enddo + ! copy in partition (gases) + do icol = 1,ncol + gicol = gicol_clr(icol + cols - 1) + colh2o(:,icol) = gh2ovmr(gicol,1:nlay) + colco2(:,icol) = gco2vmr(gicol,1:nlay) + colo3 (:,icol) = go3vmr (gicol,1:nlay) + colch4(:,icol) = gch4vmr(gicol,1:nlay) + colo2 (:,icol) = go2vmr (gicol,1:nlay) + end do + + else + + ! ------------------ + ! Cloudy gridcolumns + ! ------------------ - ! copy in partition (aerosols) - if (iaer == 10) then do icol = 1,ncol gicol = gicol_cld(icol + cols - 1) - do ibnd = 1,nbndsw - taua(1:nlay,ibnd,icol) = gtauaer(gicol,1:nlay,ibnd) - asya(1:nlay,ibnd,icol) = gasmaer(gicol,1:nlay,ibnd) - omga(1:nlay,ibnd,icol) = gssaaer(gicol,1:nlay,ibnd) + + ! assign surface albedos to bands + + ! near IR bands 14=nbndsw and 1-8 + ! 820-12850 cm-1, 0.778-12.2 um + do ibnd=1,8 + albdir(ibnd,icol) = galdir(gicol) + albdif(ibnd,icol) = galdif(gicol) + enddo + albdir(nbndsw,icol) = galdir(gicol) + albdif(nbndsw,icol) = galdif(gicol) + + ! UV/Vis bands 10-13 + ! 16000-50000 cm-1, 0.200-0.625 um + do ibnd=10,13 + albdir(ibnd,icol) = gasdir(gicol) + albdif(ibnd,icol) = gasdif(gicol) + enddo + + ! Transition band 9 + ! 12850-16000 cm-1, 0.625-0.778 um + ! Take average, dmlee + albdir(9,icol) = (gasdir(gicol)+galdir(gicol))/2. + albdif(9,icol) = (gasdif(gicol)+galdif(gicol))/2. + + enddo + + ! copy in partition (general and cloud physical props) + do icol = 1,ncol + gicol = gicol_cld(icol + cols - 1) + + play(:,icol) = gplay(gicol,1:nlay) + plev(:,icol) = gplev(gicol,1:nlay+1) + tlay(:,icol) = gtlay(gicol,1:nlay) + cld (:,icol) = gcld (gicol,1:nlay) + ciwp(:,icol) = gciwp(gicol,1:nlay) + clwp(:,icol) = gclwp(gicol,1:nlay) + rei (:,icol) = grei (gicol,1:nlay) + rel (:,icol) = grel (gicol,1:nlay) + zm (:,icol) = gzm (gicol,1:nlay) + alat (icol) = galat (gicol) + coszen(icol) = gcoszen(gicol) + enddo + + ! copy in partition (aerosols) + if (iaer == 10) then + do icol = 1,ncol + gicol = gicol_cld(icol + cols - 1) + do ibnd = 1,nbndsw + taua(1:nlay,ibnd,icol) = gtauaer(gicol,1:nlay,ibnd) + asya(1:nlay,ibnd,icol) = gasmaer(gicol,1:nlay,ibnd) + omga(1:nlay,ibnd,icol) = gssaaer(gicol,1:nlay,ibnd) + end do end do - end do - endif + endif - ! copy in partition (gases) - do icol = 1,ncol - gicol = gicol_cld(icol + cols - 1) - colh2o(:,icol) = gh2ovmr(gicol,1:nlay) - colco2(:,icol) = gco2vmr(gicol,1:nlay) - colo3 (:,icol) = go3vmr (gicol,1:nlay) - colch4(:,icol) = gch4vmr(gicol,1:nlay) - colo2 (:,icol) = go2vmr (gicol,1:nlay) - enddo + ! copy in partition (gases) + do icol = 1,ncol + gicol = gicol_cld(icol + cols - 1) + colh2o(:,icol) = gh2ovmr(gicol,1:nlay) + colco2(:,icol) = gco2vmr(gicol,1:nlay) + colo3 (:,icol) = go3vmr (gicol,1:nlay) + colch4(:,icol) = gch4vmr(gicol,1:nlay) + colo2 (:,icol) = go2vmr (gicol,1:nlay) + enddo - end if ! clear or cloudy gridcolumns + end if ! clear or cloudy gridcolumns - call MAPL_TimerOff(MAPL,"---RRTMG_PART",__RC__) + call MAPL_TimerOff(MAPL,"---RRTMG_PART",__RC__) - ! limit tiny cosine zenith angles - do icol = 1,ncol - cossza(icol) = max(zepzen,coszen(icol)) - enddo + ! limit tiny cosine zenith angles + do icol = 1,ncol + cossza(icol) = max(zepzen,coszen(icol)) + enddo - ! evaluate dry air molecules/cm^2 - ! (see details in rrtmg_lw_rad()) - do icol = 1,ncol - do ilay = 1,nlay - coldry(ilay,icol) = (plev(ilay,icol)-plev(ilay+1,icol)) * 1.e3 * avogad / & - (1.e2 * grav * ((1.-colh2o(ilay,icol)) * amd + colh2o(ilay,icol) * amw) * & - (1. + colh2o(ilay,icol))) + ! evaluate dry air molecules/cm^2 + ! (see details in rrtmg_lw_rad()) + do icol = 1,ncol + do ilay = 1,nlay + coldry(ilay,icol) = (plev(ilay,icol)-plev(ilay+1,icol)) * 1.e3 * avogad / & + (1.e2 * grav * ((1.-colh2o(ilay,icol)) * amd + colh2o(ilay,icol) * amw) * & + (1. + colh2o(ilay,icol))) + enddo enddo - enddo - ! gases also to molecules/cm^2 - do icol = 1,ncol - do ilay = 1,nlay - colh2o(ilay,icol) = coldry(ilay,icol) * colh2o(ilay,icol) - colco2(ilay,icol) = coldry(ilay,icol) * colco2(ilay,icol) - colo3 (ilay,icol) = coldry(ilay,icol) * colo3 (ilay,icol) - colch4(ilay,icol) = coldry(ilay,icol) * colch4(ilay,icol) - colo2 (ilay,icol) = coldry(ilay,icol) * colo2 (ilay,icol) + ! gases also to molecules/cm^2 + do icol = 1,ncol + do ilay = 1,nlay + colh2o(ilay,icol) = coldry(ilay,icol) * colh2o(ilay,icol) + colco2(ilay,icol) = coldry(ilay,icol) * colco2(ilay,icol) + colo3 (ilay,icol) = coldry(ilay,icol) * colo3 (ilay,icol) + colch4(ilay,icol) = coldry(ilay,icol) * colch4(ilay,icol) + colo2 (ilay,icol) = coldry(ilay,icol) * colo2 (ilay,icol) + end do end do - end do - ! cloudy gridcolumns - if (cc == 2) then - ! McICA subcolumn generation + ! We have separate loops here because MAPL profiling + ! timers must be called on all branches and all processes + ! and we do not want the timers in the if-block + call MAPL_TimerOn(MAPL,"---RRTMG_CLDSGEN",__RC__) - call generate_stochastic_clouds( & - pncol, ncol, ngptsw, nlay, & - zm, alat, dyofyr, & - play, cld, ciwp, clwp, 1.e-20, & - cldymcl, ciwpmcl, clwpmcl, & - seed_order=[4,3,2,1]) - - ! for super-layer cloud fractions - call clearCounts_threeBand( & - pncol, ncol, ngptsw, nlay, cloudLM, cloudMH, cldymcl, & - p_clearCounts) + ! cloudy gridcolumns + if (cc == 2) then + ! McICA subcolumn generation + call generate_stochastic_clouds( & + pncol, ncol, ngptsw, nlay, & + zm, alat, dyofyr, & + play, cld, ciwp, clwp, 1.e-20, & + cldymcl, ciwpmcl, clwpmcl, & + seed_order=[4,3,2,1]) + + ! for super-layer cloud fractions + call clearCounts_threeBand( & + pncol, ncol, ngptsw, nlay, cloudLM, cloudMH, cldymcl, & + p_clearCounts) + end if call MAPL_TimerOff(MAPL,"---RRTMG_CLDSGEN",__RC__) - ! cloud optical property generation call MAPL_TimerOn(MAPL,"---RRTMG_CLDPRMC",__RC__) - call cldprmc_sw( & - pncol, ncol, nlay, iceflgsw, liqflgsw, & - cldymcl, ciwpmcl, clwpmcl, rei, rel, & - taormc, taucmc, ssacmc, asmcmc) + if (cc == 2) then + ! cloud optical property generation + call cldprmc_sw( & + pncol, ncol, nlay, iceflgsw, liqflgsw, & + cldymcl, ciwpmcl, clwpmcl, rei, rel, & + taormc, taucmc, ssacmc, asmcmc) + end if call MAPL_TimerOff(MAPL,"---RRTMG_CLDPRMC",__RC__) - end if - ! Calculate information needed by the radiative transfer routine - ! that is specific to this atmosphere, especially some of the - ! coefficients and indices needed to compute the optical depths - ! by interpolating data from stored reference atmospheres. + ! Calculate information needed by the radiative transfer routine + ! that is specific to this atmosphere, especially some of the + ! coefficients and indices needed to compute the optical depths + ! by interpolating data from stored reference atmospheres. + + call MAPL_TimerOn(MAPL,"---RRTMG_SETCOEF",__RC__) + call setcoef_sw( & + pncol, ncol, nlay, play, tlay, coldry, & + colch4, colco2, colh2o, colmol, colo2, colo3, & + laytrop, jp, jt, jt1, fac00, fac01, fac10, fac11, & + selffac, selffrac, indself, forfac, forfrac, indfor) + call MAPL_TimerOff(MAPL,"---RRTMG_SETCOEF",__RC__) + + ! compute sw radiative fluxes + call spcvmc_sw(MAPL, & + cc, pncol, ncol, nlay, & + albdif, albdir, & + cldymcl, taucmc, asmcmc, ssacmc, taormc, & + taua, asya, omga, cossza, adjflux, & + isolvar, svar_f, svar_s, svar_i, & + svar_f_bnd, svar_s_bnd, svar_i_bnd, & + laytrop, jp, jt, jt1, & + colch4, colco2, colh2o, colmol, colo2, colo3, & + fac00, fac01, fac10, fac11, & + cloudLM, cloudMH, & + selffac, selffrac, indself, forfac, forfrac, indfor, & + zbbfd, zbbfu, zbbcd, zbbcu, zuvfd, zuvcd, znifd, znicd, & + zbbfddir, zbbcddir, zuvfddir, zuvcddir, znifddir, znicddir,& + znirr, znirf, zparr, zparf, zuvrr, zuvrf, fndsbnd, & + ztautp, ztauhp, ztaump, ztaulp, & + do_drfband, zdrband, zdfband, & + __RC__) + + ! Copy out up and down, clear- and all-sky fluxes to output arrays. + ! Vertical indexing goes from bottom to top; reverse here for GCM if necessary. + + call MAPL_TimerOn(MAPL,"---RRTMG_PART",__RC__) + + if (cc == 1) then ! clear gridcolumns - call MAPL_TimerOn(MAPL,"---RRTMG_SETCOEF",__RC__) - call setcoef_sw( & - pncol, ncol, nlay, play, tlay, coldry, & - colch4, colco2, colh2o, colmol, colo2, colo3, & - laytrop, jp, jt, jt1, fac00, fac01, fac10, fac11, & - selffac, selffrac, indself, forfac, forfrac, indfor) - call MAPL_TimerOff(MAPL,"---RRTMG_SETCOEF",__RC__) + do icol = 1,ncol + gicol = gicol_clr(icol + cols - 1) - ! compute sw radiative fluxes - call spcvmc_sw(MAPL, & - cc, pncol, ncol, nlay, & - albdif, albdir, & - cldymcl, taucmc, asmcmc, ssacmc, taormc, & - taua, asya, omga, cossza, adjflux, & - isolvar, svar_f, svar_s, svar_i, & - svar_f_bnd, svar_s_bnd, svar_i_bnd, & - laytrop, jp, jt, jt1, & - colch4, colco2, colh2o, colmol, colo2, colo3, & - fac00, fac01, fac10, fac11, & - cloudLM, cloudMH, & - selffac, selffrac, indself, forfac, forfrac, indfor, & - zbbfd, zbbfu, zbbcd, zbbcu, zuvfd, zuvcd, znifd, znicd, & - zbbfddir, zbbcddir, zuvfddir, zuvcddir, znifddir, znicddir,& - znirr, znirf, zparr, zparf, zuvrr, zuvrf, fndsbnd, & - ztautp, ztauhp, ztaump, ztaulp, & - do_drfband, zdrband, zdfband, & - __RC__) - - ! Copy out up and down, clear- and all-sky fluxes to output arrays. - ! Vertical indexing goes from bottom to top; reverse here for GCM if necessary. + ! super-layer clear counts + do n = 1,4 + clearCounts (gicol,n) = ngptsw + end do - call MAPL_TimerOn(MAPL,"---RRTMG_PART",__RC__) + ! up and down fluxes + do ilev = 1,nlay+1 + swuflxc(gicol,ilev) = zbbcu(ilev,icol) + swdflxc(gicol,ilev) = zbbcd(ilev,icol) + swuflx (gicol,ilev) = zbbfu(ilev,icol) + swdflx (gicol,ilev) = zbbfd(ilev,icol) + enddo - if (cc == 1) then ! clear gridcolumns + ! super-layer optical thicknesses + tautp(gicol) = 0. + tauhp(gicol) = 0. + taump(gicol) = 0. + taulp(gicol) = 0. - do icol = 1,ncol - gicol = gicol_clr(icol + cols - 1) - - ! super-layer clear counts - do n = 1,4 - clearCounts (gicol,n) = ngptsw - end do - - ! up and down fluxes - do ilev = 1,nlay+1 - swuflxc(gicol,ilev) = zbbcu(ilev,icol) - swdflxc(gicol,ilev) = zbbcd(ilev,icol) - swuflx (gicol,ilev) = zbbfu(ilev,icol) - swdflx (gicol,ilev) = zbbfd(ilev,icol) enddo - ! super-layer optical thicknesses - tautp(gicol) = 0. - tauhp(gicol) = 0. - taump(gicol) = 0. - taulp(gicol) = 0. - - enddo - - ! surface broadband fluxes - do icol = 1,ncol - gicol = gicol_clr(icol + cols - 1) - nirr(gicol) = znirr(icol) - nirf(gicol) = znirf(icol) - znirr(icol) - parr(gicol) = zparr(icol) - parf(gicol) = zparf(icol) - zparr(icol) - uvrr(gicol) = zuvrr(icol) - uvrf(gicol) = zuvrf(icol) - zuvrr(icol) - end do - - ! net downward flux at surface in bands - do ibnd = 1,nbndsw + ! surface broadband fluxes do icol = 1,ncol gicol = gicol_clr(icol + cols - 1) - fswband(gicol,ibnd) = fndsbnd(icol,ibnd) + nirr(gicol) = znirr(icol) + nirf(gicol) = znirf(icol) - znirr(icol) + parr(gicol) = zparr(icol) + parf(gicol) = zparf(icol) - zparr(icol) + uvrr(gicol) = zuvrr(icol) + uvrf(gicol) = zuvrf(icol) - zuvrr(icol) end do - end do - ! downward beam and diffuse fluxes at surface in bands - if (do_drfband) then + ! net downward flux at surface in bands do ibnd = 1,nbndsw do icol = 1,ncol gicol = gicol_clr(icol + cols - 1) - drband(gicol,ibnd) = zdrband(icol,ibnd) - dfband(gicol,ibnd) = zdfband(icol,ibnd) + fswband(gicol,ibnd) = fndsbnd(icol,ibnd) end do end do - end if - else ! cloudy columns + ! downward beam and diffuse fluxes at surface in bands + if (do_drfband) then + do ibnd = 1,nbndsw + do icol = 1,ncol + gicol = gicol_clr(icol + cols - 1) + drband(gicol,ibnd) = zdrband(icol,ibnd) + dfband(gicol,ibnd) = zdfband(icol,ibnd) + end do + end do + end if - do icol = 1,ncol - gicol = gicol_cld(icol + cols - 1) - do n = 1,4 - clearCounts (gicol,n) = p_clearCounts(n,icol) - end do - do ilev = 1,nlay+1 - swuflxc(gicol,ilev) = zbbcu(ilev,icol) - swdflxc(gicol,ilev) = zbbcd(ilev,icol) - swuflx (gicol,ilev) = zbbfu(ilev,icol) - swdflx (gicol,ilev) = zbbfd(ilev,icol) - enddo - tautp(gicol) = ztautp(icol) - tauhp(gicol) = ztauhp(icol) - taump(gicol) = ztaump(icol) - taulp(gicol) = ztaulp(icol) - enddo + else ! cloudy columns - do icol = 1,ncol - gicol = gicol_cld(icol + cols - 1) - nirr(gicol) = znirr(icol) - nirf(gicol) = znirf(icol) - znirr(icol) - parr(gicol) = zparr(icol) - parf(gicol) = zparf(icol) - zparr(icol) - uvrr(gicol) = zuvrr(icol) - uvrf(gicol) = zuvrf(icol) - zuvrr(icol) - enddo + do icol = 1,ncol + gicol = gicol_cld(icol + cols - 1) + do n = 1,4 + clearCounts (gicol,n) = p_clearCounts(n,icol) + end do + do ilev = 1,nlay+1 + swuflxc(gicol,ilev) = zbbcu(ilev,icol) + swdflxc(gicol,ilev) = zbbcd(ilev,icol) + swuflx (gicol,ilev) = zbbfu(ilev,icol) + swdflx (gicol,ilev) = zbbfd(ilev,icol) + enddo + tautp(gicol) = ztautp(icol) + tauhp(gicol) = ztauhp(icol) + taump(gicol) = ztaump(icol) + taulp(gicol) = ztaulp(icol) + enddo - ! net downward flux at surface in bands - do ibnd = 1,nbndsw do icol = 1,ncol gicol = gicol_cld(icol + cols - 1) - fswband(gicol,ibnd) = fndsbnd(icol,ibnd) - end do - end do + nirr(gicol) = znirr(icol) + nirf(gicol) = znirf(icol) - znirr(icol) + parr(gicol) = zparr(icol) + parf(gicol) = zparf(icol) - zparr(icol) + uvrr(gicol) = zuvrr(icol) + uvrf(gicol) = zuvrf(icol) - zuvrr(icol) + enddo - ! downward beam and diffuse fluxes at surface in bands - if (do_drfband) then + ! net downward flux at surface in bands do ibnd = 1,nbndsw do icol = 1,ncol gicol = gicol_cld(icol + cols - 1) - drband(gicol,ibnd) = zdrband(icol,ibnd) - dfband(gicol,ibnd) = zdfband(icol,ibnd) + fswband(gicol,ibnd) = fndsbnd(icol,ibnd) end do end do - end if - endif ! clear/cloudy + ! downward beam and diffuse fluxes at surface in bands + if (do_drfband) then + do ibnd = 1,nbndsw + do icol = 1,ncol + gicol = gicol_cld(icol + cols - 1) + drband(gicol,ibnd) = zdrband(icol,ibnd) + dfband(gicol,ibnd) = zdfband(icol,ibnd) + end do + end do + end if + + endif ! clear/cloudy - call MAPL_TimerOff(MAPL,"---RRTMG_PART",__RC__) + call MAPL_TimerOff(MAPL,"---RRTMG_PART",__RC__) - enddo ! over partitions + enddo ! over partitions + + endif ! npart > 0 enddo ! outer loop (cc) over clear then cloudy columns @@ -1305,7 +1340,7 @@ real function earth_sun(idn) ! ! Purpose: Function to calculate the correction factor of Earth's orbit ! for current day of the year - ! + ! ! idn : Day of the year ! earth_sun : square of the ratio of mean to actual Earth-Sun distance !----------------------------------------------------------------------- @@ -1316,7 +1351,7 @@ real function earth_sun(idn) real :: gamma - gamma = 2. * pi * (idn-1)/365. + gamma = 2. * pi * (idn-1)/365. ! Use Iqbal's equation 1.2.1