Changeset 2559

Timestamp:

Sep 7, 2021, 12:10:35 PM (3 years ago)

Author:

emillour

Message:

Mars GCM:
Some code cleanup and refactoring around wstats:

• turn wstats.F90 in a module
• remove no useless statto_mod.F90
• incorporate auxilliary inistats and mkstats routines in wstats_mod.F90
• move flag "callstats" from callkeys.h to being a module variable of wstats_mod

EM

Location:

trunk/LMDZ.MARS

Files:

• README (modified) (3 diffs)
• libf/aeronomars/calchim_mod.F90 (modified) (1 diff)
• libf/aeronomars/surfacearea.F (modified) (1 diff)
• libf/phymars/callkeys.h (modified) (2 diffs)
• libf/phymars/conf_phys.F (modified) (1 diff)
• libf/phymars/inistats.F (deleted)
• libf/phymars/mkstat.F90 (deleted)
• libf/phymars/physiq_mod.F (modified) (1 diff)
• libf/phymars/statto_mod.F90 (deleted)
• libf/phymars/wstats_mod.F90 (moved) (moved from trunk/LMDZ.MARS/libf/phymars/wstats.F90) (3 diffs)

trunk/LMDZ.MARS/README

@@ -2466 +2466 @@
 microphysics (just like water clouds). The amplitude is spantCO2, also read from callphys.def
 - The previous logical microphysco2 has been removed.
-- Cloud opacity at 1µm is now computed in the co2cloud.F routine
+- Cloud opacity at 1µm is now computed in the co2cloud.F routine
 - Most of the co2 ice clouds scheme writediagfi are now in co2clouds.F
 - Some cleaning of the CO2 ice clouds routine has been done. Not perfect yet!
@@ -2482 +2482 @@
   if a gravity wave propagates to the mesosphere or not (wheter the wave saturates on its way or not). 
   The sub-grid temperature distribution is cancelled if the saturation index value of the column (between 12 and 80 km) is > 0.1, and is applied (+/- 3K) otherwise. if the keyword satindexco2=.true. in callphys.def (only applies if CLFvaryingCO2=.true. anyway). If set to .false., deactivate this filter and the sub grid T distribution applies everywhere and everytime. See comments in co2clouds.F
-- All that you need to launch a GCM run with co2 clouds is described in co2clouds.F comments. Ehouarn Millour and Anni Määttänen have the files. 
+- All that you need to launch a GCM run with co2 clouds is described in co2clouds.F comments. Ehouarn Millour and Anni Määttänen have the files. 
  
 == 15/11/2017 == EM
@@ -3443 +3443 @@
 (between levels with no missing values in both files, which vary with space and time) and outputs
 the ratio tau_GCM/tau_MCS to be used to normalize the GCM dust opacity profiles when comparing to MCS
+
+== 07/09/2021 == EM
+Some code cleanup and refactoring around wstats:
+- turn wstats.F90 in a module
+- remove no useless statto_mod.F90
+- incorporate auxilliary inistats and mkstats routines in wstats_mod.F90
+- move flag "callstats" from callkeys.h to being a module variable of wstats_mod

trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90

@@ -31 +31 @@
       use photolysis_mod, only: init_photolysis, nphot
       use iono_h, only: temp_elect
+      use wstats_mod, only: callstats, wstats
 
       implicit none

trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F

@@ -8 +8 @@
      &                      igcm_ccn_number, varian, ccn_factor
       use conc_mod, only: rnew
-      USE comcstfi_h
+      use comcstfi_h, only: pi
+      use wstats_mod, only: callstats, wstats
       implicit none
 

trunk/LMDZ.MARS/libf/phymars/callkeys.h

@@ -7 +7 @@
 
       COMMON/callkeys_l/callrad,calldifv,calladj,callcond,callsoil      &
-     &   ,season,diurnal,lwrite,calllott,callstats,calleofdump          &
+     &   ,season,diurnal,lwrite,calllott,calleofdump                    &
      &   ,callnirco2,callnlte,callthermos,callconduct,calleuv           &
      &   ,callmolvis,callmoldiff,thermochem,thermoswater,callemis       &
@@ -28 +28 @@
       LOGICAL callrad,calldifv,calladj,callcond,callsoil,               &
      &   season,diurnal,lwrite,calllott,calllott_nonoro                 &
-     &   ,callstats,calleofdump                                         &
+     &   ,calleofdump                                                   &
      &   ,callnirco2,callnlte,callthermos,callconduct,                  &
      &    calleuv,callmolvis,callmoldiff,thermochem,thermoswater        &

trunk/LMDZ.MARS/libf/phymars/conf_phys.F

@@ -43 +43 @@
      &                      ini_scatterers,tauvis
       use datafile_mod, only: datadir
+      use wstats_mod, only: callstats
       use calchim_mod, only: ichemistry
       use co2condens_mod, only: scavco2cond

trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

@@ -84 +84 @@
       use phyredem, only: physdem0, physdem1
       use phyetat0_mod, only: phyetat0, tab_cntrl_mod
+      use wstats_mod, only: callstats, wstats, mkstats
       use eofdump_mod, only: eofdump
       USE vertical_layers_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt

trunk/LMDZ.MARS/libf/phymars/wstats_mod.F90

@@ +1 @@
+module wstats_mod
+
+! module containing parameters and routines to generate the "stats.nc" file
+! which will contain a mean statistical day of variables extracted at set
+! times of day every day of the simulation, and the standard deviations thereof
+
+implicit none
+
+logical,save :: callstats ! global flag to trigger generating a stats.nc
+                          ! file or not. Initialized in conf_gcm()
+!$OMP THREADPRIVATE(callstats)
+
+integer,save :: istats ! calculate stats every istats physics timestep,
+                       ! starting at first call. Initialized by inistats()
+!$OMP THREADPRIVATE(istats)
+
+integer,parameter :: istime=12 ! number of time steps per sol to store
+
+integer,save :: count(istime) ! count tab to know the variable record
+!$OMP THREADPRIVATE(istats,count)
+
+contains
+
 subroutine wstats(ngrid,nom,titre,unite,dim,px)
 
-use statto_mod, only: istats,istime,count
+! Main routine to call to trigger writing a given field to the "stats.nc" file 
+
 use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather, klon_mpi_begin
 use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
@@ -335 +359 @@
 end subroutine wstats
 
-!======================================================
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine inistats(ierr)
+
+! routine to initialize the stats file (i.e. create file, write
+! time independent coordinates, etc.)
+
+use mod_phys_lmdz_para, only : is_master
+USE vertical_layers_mod, ONLY: ap,bp,aps,bps,preff,pseudoalt,presnivs
+USE nrtype, ONLY: pi
+USE time_phylmdz_mod, ONLY: daysec,dtphys
+USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
+implicit none
+
+include "netcdf.inc"
+
+integer,intent(out) :: ierr
+
+integer :: nid
+integer :: l,nsteppd
+real, dimension(nbp_lev) ::  sig_s
+real,allocatable :: lon_reg_ext(:) ! extended longitudes
+integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time
+real, dimension(istime) :: lt
+integer :: nvarid
+
+
+IF (nbp_lon*nbp_lat==1) THEN
+  ! 1D model
+  ALLOCATE(lon_reg_ext(1))
+ELSE
+  ! 3D model
+  ALLOCATE(lon_reg_ext(nbp_lon+1))
+ENDIF
+
+write (*,*) 
+write (*,*) '                        || STATS ||'
+write (*,*) 
+write (*,*) 'daysec',daysec
+write (*,*) 'dtphys',dtphys
+nsteppd=nint(daysec/dtphys)
+write (*,*) 'nsteppd=',nsteppd
+
+if (abs(float(nsteppd)-daysec/dtphys).gt.1.e-8*daysec) then
+  call abort_physic("inistats",'1 sol .ne. n physics steps',1)
+endif
+
+if(mod(nsteppd,istime).ne.0) then
+  call abort_physic("inistats",'1 sol .ne. n*istime physics steps',1)
+endif
+
+istats=nsteppd/istime
+write (*,*) 'istats=',istats
+write (*,*) 'Storing ',istime,'times per day'
+write (*,*) 'thus every ',istats,'physical timestep '
+write (*,*) 
+
+do l= 1, nbp_lev
+  sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
+  pseudoalt(l)=-10.*log(presnivs(l)/preff)   
+enddo
+
+lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
+IF (nbp_lon*nbp_lat/=1) THEN
+  ! In 3D, add extra redundant point (180 degrees,
+  ! since lon_reg starts at -180)
+  lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+ENDIF
+
+if (is_master) then
+  ! only the master needs do this
+  ierr = NF_CREATE("stats.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
+  if (ierr.ne.NF_NOERR) then
+    write (*,*) NF_STRERROR(ierr)
+    call abort_physic("inistats",'failed creating stats.nc',1)
+  endif
+
+  ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_lat)
+  IF (nbp_lon*nbp_lat==1) THEN
+    ierr = NF_DEF_DIM (nid, "longitude", 1, idim_lon)
+  ELSE
+    ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_lon)
+  ENDIF
+  ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
+  ierr = NF_DEF_DIM (nid, "llmp1", nbp_lev+1, idim_llmp1)
+  ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_time)
+
+  ierr = NF_ENDDEF(nid)
+
+  call def_var_stats(nid,"Time","Time", &
+                     "days since 0000-00-0 00:00:00",1, &
+                     [idim_time],nvarid,ierr)
+  ! Time is initialised later by mkstats subroutine
+
+  call def_var_stats(nid,"latitude","latitude", &
+                     "degrees_north",1,[idim_lat],nvarid,ierr)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
+#else
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
+#endif
+
+  call def_var_stats(nid,"longitude","East longitude", &
+                     "degrees_east",1,[idim_lon],nvarid,ierr)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
+#else
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
+#endif
+
+  ! Niveaux verticaux, aps et bps
+  ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+  ierr = NF_DEF_VAR (nid,"altitude", NF_DOUBLE, 1,idim_llm,nvarid)
+#else
+  ierr = NF_DEF_VAR (nid,"altitude", NF_FLOAT, 1,idim_llm,nvarid)
+#endif
+  ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",8,"altitude")
+  ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km")
+  ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up")
+  ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt)
+#else
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
+#endif 
+  call def_var_stats(nid,"aps","hybrid pressure at midlayers", &
+                     " ",1,[idim_llm],nvarid,ierr)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
+#else
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,aps)
+#endif
+
+  call def_var_stats(nid,"bps","hybrid sigma at midlayers", &
+                     " ",1,[idim_llm],nvarid,ierr)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
+#else 
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,bps)
+#endif
+
+  ierr=NF_CLOSE(nid)
+
+endif ! of if (is_master)
+
+end subroutine inistats
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
 subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
+
+! routine to initialize writing a variable in the stats file
+
 use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
 
@@ -516 +693 @@
 end subroutine def_var_stats
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine mkstats(ierr)
+
+!  This routine finalizes tha stats.nc file from sums and sums of squares
+!  to means and standard deviations. The data file is overwritten in place.
+
+use mod_phys_lmdz_para, only : is_master
+use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
+
+implicit none
+
+include "netcdf.inc"
+
+integer,intent(out) :: ierr ! netCDF return error code
+integer :: nid,nbvar,i,ndims,lt,nvarid
+integer, dimension(4) :: id,varid,start,size
+integer, dimension(5) :: dimids
+character (len=50) :: name,nameout,units,title
+real,allocatable :: sum3d(:,:,:),square3d(:,:,:),mean3d(:,:,:),sd3d(:,:,:)
+real,allocatable :: sum2d(:,:),square2d(:,:),mean2d(:,:),sd2d(:,:)
+real, dimension(istime) :: time
+real, dimension(nbp_lat) :: lat
+real,allocatable :: lon(:)
+real, dimension(nbp_lev) :: alt
+logical :: lcopy=.true.
+!integer :: latid,lonid,altid,timeid
+integer :: meanid,sdid
+!integer, dimension(4) :: dimout
+
+if (is_master) then
+! only the master needs do all this
+
+! Incrementation of count for the last step, which is not done in wstats
+count(istime)=count(istime)+1
+
+if (klon_glo>1) then
+  allocate(lon(nbp_lon+1))
+  allocate(sum3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(square3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(mean3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(sd3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(sum2d(nbp_lon+1,nbp_lat))
+  allocate(square2d(nbp_lon+1,nbp_lat))
+  allocate(mean2d(nbp_lon+1,nbp_lat))
+  allocate(sd2d(nbp_lon+1,nbp_lat))
+else ! 1D model case
+  allocate(lon(1))
+  allocate(sum3d(1,1,nbp_lev))
+  allocate(square3d(1,1,nbp_lev))
+  allocate(mean3d(1,1,nbp_lev))
+  allocate(sd3d(1,1,nbp_lev))
+  allocate(sum2d(1,1))
+  allocate(square2d(1,1))
+  allocate(mean2d(1,1))
+  allocate(sd2d(1,1))
+endif
+
+ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
+
+! We catch the id of dimensions of the stats file
+
+ierr= NF_INQ_DIMID(nid,"latitude",id(1))
+ierr= NF_INQ_DIMID(nid,"longitude",id(2))
+ierr= NF_INQ_DIMID(nid,"altitude",id(3))
+ierr= NF_INQ_DIMID(nid,"Time",id(4))
+
+ierr= NF_INQ_VARID(nid,"latitude",varid(1))
+ierr= NF_INQ_VARID(nid,"longitude",varid(2))
+ierr= NF_INQ_VARID(nid,"altitude",varid(3))
+ierr= NF_INQ_VARID(nid,"Time",varid(4))
+
+! Time initialisation
+
+do i=1,istime
+   time(i)=i*24./istime
+#ifdef NC_DOUBLE
+   ierr= NF_PUT_VARA_DOUBLE(nid,varid(4),i,1,time(i))
+#else
+   ierr= NF_PUT_VARA_REAL(nid,varid(4),i,1,time(i))
+#endif
+enddo
+
+! We catche the values of the variables
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,varid(1),lat)
+         ierr = NF_GET_VAR_DOUBLE(nid,varid(2),lon)
+         ierr = NF_GET_VAR_DOUBLE(nid,varid(3),alt)
+#else
+         ierr = NF_GET_VAR_REAL(nid,varid(1),lat)
+         ierr = NF_GET_VAR_REAL(nid,varid(2),lon)
+         ierr = NF_GET_VAR_REAL(nid,varid(3),alt)
+#endif
+
+! We catch the number of variables in the stats file
+ierr = NF_INQ_NVARS(nid,nbvar)
+
+! to catche the "real" number of variables (without the "additionnal variables")
+nbvar=(nbvar-4)/2 
+
+do i=1,nbvar
+   varid=(i-1)*2+5
+
+   ! What's the variable's name?
+   ierr=NF_INQ_VARNAME(nid,varid,name)
+   write(*,*) "OK variable ",name
+   ! Its units?
+   units=" "
+   ierr=NF_GET_ATT_TEXT(nid,varid,"units",units)
+   ! Its title?
+   title=" "
+   ierr=NF_GET_ATT_TEXT(nid,varid,"title",title)
+   ! Its number of dimensions?   
+   ierr=NF_INQ_VARNDIMS(nid,varid,ndims)
+   ! Its values?
+
+   if(ndims==4) then ! lat, lon, alt & time
+
+!      dimout(1)=lonid
+!      dimout(2)=latid
+!      dimout(3)=altid
+!      dimout(4)=timeid
+
+      if (klon_glo>1) then ! general case
+        size=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+      else ! 1D model
+        size=(/1,1,nbp_lev,1/)
+      endif
+      do lt=1,istime
+         start=(/1,1,1,lt/)
+         ! Extraction of the "source" variables
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(nid,varid,start,size,sum3d)
+         ierr = NF_GET_VARA_DOUBLE(nid,varid+1,start,size,square3d)
+#else
+         ierr = NF_GET_VARA_REAL(nid,varid,start,size,sum3d)
+         ierr = NF_GET_VARA_REAL(nid,varid+1,start,size,square3d)
+#endif
+         ! Calculation of these variables
+         mean3d=sum3d/count(lt)
+         sd3d=sqrt(max(0.,square3d/count(lt)-mean3d**2))
+         ! Writing of the variables
+#ifdef NC_DOUBLE
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,mean3d)
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid+1,start,size,sd3d)
+#else
+         ierr = NF_PUT_VARA_REAL(nid,varid,start,size,mean3d)
+         ierr = NF_PUT_VARA_REAL(nid,varid+1,start,size,sd3d)
+#endif
+      enddo
+
+    else if (ndims.eq.3) then
+
+!      dimout(1)=lonid
+!      dimout(2)=latid
+!      dimout(3)=timeid
+
+      if (klon_glo>1) then ! general case
+        size=(/nbp_lon+1,nbp_lat,1,0/)
+      else
+        size=(/1,1,1,0/)
+      endif
+      do lt=1,istime
+         start=(/1,1,lt,0/)
+         ! Extraction of the "source" variables
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(nid,varid,start,size,sum2d)
+         ierr = NF_GET_VARA_DOUBLE(nid,varid+1,start,size,square2d)
+#else
+         ierr = NF_GET_VARA_REAL(nid,varid,start,size,sum2d)
+         ierr = NF_GET_VARA_REAL(nid,varid+1,start,size,square2d)
+#endif
+         ! Calculation of these variables
+         mean2d=sum2d/count(lt)
+         sd2d=sqrt(max(0.,square2d/count(lt)-mean2d**2))
+         ! Writing of the variables
+#ifdef NC_DOUBLE
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,mean2d)
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid+1,start,size,sd2d)
+#else
+         ierr = NF_PUT_VARA_REAL(nid,varid,start,size,mean2d)
+         ierr = NF_PUT_VARA_REAL(nid,varid+1,start,size,sd2d)
+#endif
+      enddo
+
+    endif 
+enddo
+
+ierr= NF_CLOSE(nid)
+
+endif ! of if (is_master)
+
+end subroutine mkstats
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+end module wstats_mod