modtimestat: add radiative fluxes at surface, TOM, TOA for RCEMIP

dalesteam · Nov 24, 2023 · a372524 · a372524
1 parent 7ab0999
commit a372524
Showing 1 changed file with 65 additions and 5 deletions.
diff --git a/src/modtimestat.f90 b/src/modtimestat.f90
@@ -44,6 +44,7 @@ module modtimestat
 
   integer :: nvar
   integer :: ncid,nrec = 0
+  integer :: ivar_rad = 0  ! starting index for radiation quantities
   character(80) :: fname = 'tmser.xxx.nc'
   !character(80),dimension(nvar,4) :: ncname
   character(80), allocatable, dimension(:,:)    :: ncname
@@ -82,12 +83,13 @@ module modtimestat
 contains
 !> Initializing Timestat. Read out the namelist, initializing the variables
   subroutine inittimestat
-    use modmpi,    only : myid,comm3d,mpi_logical,mpierr,mpi_integer, D_MPI_BCAST
+    use modmpi,    only : myid,comm3d,mpierr,D_MPI_BCAST
     use modglobal, only : ifnamopt, fname_options,cexpnr,dtmax,ifoutput,dtav_glob,tres,&
                           ladaptive,k1,kmax,rd,rv,dt_lim,btime,i1,j1,lwarmstart,checknamelisterror
     use modfields, only : thlprof,qtprof,svprof
     use modsurfdata, only : isurf, lhetero, xpatches, ypatches
     use modstat_nc, only : lnetcdf, open_nc, define_nc, ncinfo, nctiminfo
+    use modraddata, only : iradiation
     implicit none
     integer :: ierr,k,location = 1
     real :: gradient = 0.0
@@ -227,11 +229,15 @@ subroutine inittimestat
           endif
        endif
 
-      if (lnetcdf) then
+       if (lnetcdf) then
+        nvar = 23
+        ivar_rad = 24
         if(isurf == 1) then
-          nvar = 34
-        else
-          nvar = 23
+           nvar = nvar + 11
+           ivar_rad = ivar_rad + 11
+        end if
+        if (iradiation /= 0) then
+           nvar = nvar + 19
         end if
 
         allocate(ncname(nvar,4))
@@ -274,6 +280,31 @@ subroutine inittimestat
           call ncinfo(ncname(33,:),'rssoil','Soil evaporation resistance','s/m','time')
           call ncinfo(ncname(34,:),'rsveg','Vegitation resistance','s/m','time')
         end if
+
+        if (iradiation /= 0) then
+           call ncinfo(ncname(ivar_rad+ 0,:),'rlds',   'surface downwelling longwave flux','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+ 1,:),'rlus',   'surface upwelling longwave flux','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+ 2,:),'rsds',   'surface downwelling shortwave flux','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+ 3,:),'rsus',   'surface upwelling shortwave flux','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+ 4,:),'rsdscs', 'surface downwelling shortwave flux - clear sky','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+ 5,:),'rsuscs', 'surface upwelling shortwave flux - clear sky','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+ 6,:),'rldscs', 'surface downwelling longwave flux - clear sky','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+ 7,:),'rluscs', 'surface upwelling longwave flux - clear sky','W/m^2','time')
+
+           call ncinfo(ncname(ivar_rad+ 8,:),'rsdt',   'TOA incoming shortwave flux','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+ 9,:),'rsut',   'TOA outgoing shortwave flux','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+10,:),'rlut',   'TOA outgoing longwave flux','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+11,:),'rsutcs', 'TOA outgoing shortwave flux -clear sky','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+12,:),'rlutcs', 'TOA outgoing longwave flux -clear sky','W/m^2','time')
+
+           call ncinfo(ncname(ivar_rad+13,:),'rsdtm',  'TOM incoming shortwave flux','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+14,:),'rldtm',  'TOM incoming longwave flux','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+15,:),'rsutm',  'TOM outgoing shortwave flux','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+16,:),'rlutm',  'TOM outgoing longwave flux','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+17,:),'rsutmcs','TOM outgoing shortwave flux -clear sky','W/m^2','time')
+           call ncinfo(ncname(ivar_rad+18,:),'rlutmcs','TOM outgoing longwave flux -clear sky','W/m^2','time')
+        end if
+
         call open_nc(fname,  ncid,nrec)
         if(nrec==0) call define_nc( ncid, NVar, ncname)
       end if
@@ -350,6 +381,8 @@ subroutine timestat
     use modsurface, only : patchxnr,patchynr
     use modmpi,     only : mpi_sum,mpi_max,mpi_min,comm3d,mpierr,myid, D_MPI_ALLREDUCE
     use modstat_nc,  only : lnetcdf, writestat_nc,nc_fillvalue
+    use modraddata, only :  lwd,lwu,swd,swu,lwdca,lwuca,swdca,swuca,SW_up_TOA,SW_dn_TOA,LW_up_TOA,&
+                            SW_up_ca_TOA,LW_up_ca_TOA,iradiation
     implicit none
 
     real   :: zbaseavl, ztopavl, ztopmaxl, ztop,zbaseminl
@@ -884,6 +917,33 @@ subroutine timestat
           vars(33) = rssoilav
           vars(34) = rsvegav
         end if
+        if (iradiation /= 0) then
+           vars(ivar_rad+ 0) = abs(sum(lwd(2:i1,2:j1,1)))  !'rlds',   'surface downwelling longwave flux'
+           vars(ivar_rad+ 1) = abs(sum(lwu(2:i1,2:j1,1)))  !'rlus',   'surface upwelling longwave flux'
+           vars(ivar_rad+ 2) = abs(sum(swd(2:i1,2:j1,1)))  !'rsds',   'surface downwelling shortwave flux'
+           vars(ivar_rad+ 3) = abs(sum(swu(2:i1,2:j1,1)))  !'rsus',   'surface upwelling shortwave flux'
+
+           vars(ivar_rad+ 4) = abs(sum(swdca(2:i1,2:j1,1)))  !'rsdscs', 'surface downwelling shortwave flux - clear sky'
+           vars(ivar_rad+ 5) = abs(sum(swuca(2:i1,2:j1,1)))  !'rsuscs', 'surface upwelling shortwave flux - clear sky'
+           vars(ivar_rad+ 6) = abs(sum(lwdca(2:i1,2:j1,1)))  !'rldscs', 'surface downwelling longwave flux - clear sky'
+           vars(ivar_rad+ 7) = abs(sum(lwuca(2:i1,2:j1,1)))  !'rluscs', 'surface upwelling longwave flux - clear sky'
+
+           vars(ivar_rad+ 8) = abs(sum(SW_dn_TOA(2:i1,2:j1)))    !'rsdt',   'TOA incoming shortwave flux','W/m^2'
+           vars(ivar_rad+ 9) = abs(sum(SW_up_TOA(2:i1,2:j1)))    !'rsut',   'TOA outgoing shortwave flux','W/m^2'
+           vars(ivar_rad+10) = abs(sum(LW_up_TOA(2:i1,2:j1)))    !'rlut',   'TOA outgoing longwave flux','W/m^2'
+           vars(ivar_rad+11) = abs(sum(SW_up_ca_TOA(2:i1,2:j1))) !'rsutcs', 'TOA outgoing shortwave flux -clear sky'
+           vars(ivar_rad+12) = abs(sum(LW_up_ca_TOA(2:i1,2:j1))) !'rlutcs', 'TOA outgoing longwave flux -clear sky'
+
+           vars(ivar_rad+13) = abs(sum(swd(2:i1,2:j1,kmax)))   !'rsdtm',  'TOM incoming shortwave flux'
+           vars(ivar_rad+14) = abs(sum(lwd(2:i1,2:j1,kmax)))   !'rsdtm',  'TOM incoming longwave flux'
+           vars(ivar_rad+15) = abs(sum(swu(2:i1,2:j1,kmax)))   !'rsutm',  'TOM outgoing shortwave flux'
+           vars(ivar_rad+16) = abs(sum(lwu(2:i1,2:j1,kmax)))   !'rlutm',  'TOM outgoing longwave flux'
+           vars(ivar_rad+17) = abs(sum(swuca(2:i1,2:j1,kmax))) !'rsutmcs','TOM outgoing shortwave flux -clear sky'
+           vars(ivar_rad+18) = abs(sum(lwuca(2:i1,2:j1,kmax))) !'rlutmcs','TOM outgoing longwave flux -clear sky'
+
+           call D_MPI_ALLREDUCE(vars(ivar_rad:ivar_rad+18), 19,  MPI_SUM, comm3d,mpierr)
+           vars(ivar_rad:ivar_rad+18) = vars(ivar_rad:ivar_rad+18) / ijtot
+        end if
 
         call writestat_nc(ncid,nvar,ncname,vars,nrec,.true.)
       end if