Changeset 2238 for trunk/LMDZ.MARS/util

Timestamp:

Feb 4, 2020, 6:46:34 PM (5 years ago)

Author:

abierjon

Message:

Mars GCM:
Small improvements in simu_MCS_temp.F90 : 1) The binning method is now more consistent with the MCS file ; 2) stats files are more rigorously recognized.
AB

File:

• trunk/LMDZ.MARS/util/simu_MCS_temp.F90 (modified) (27 diffs)

trunk/LMDZ.MARS/util/simu_MCS_temp.F90

@@ -8 +8 @@
 !===================================================================================================
 ! This program loads Observer's data file (MCS-type, binned by Luca Montabone), serving as a 
-! reference, then extract the dayside and nightside temperature data from either :
+! reference, then extracts the dayside and nightside temperature data from either :
 !       - a GCM simulation (diagfi.nc) that covers the observation period,
 !       - a GCM stats file (stats.nc) after extending it to cover the observation period.
@@ -15 +15 @@
 ! the Observer's spatial coordinates, at every GCM sol value in each 5°-interval of Ls, before 
 ! being averaged to create output bins for each interval. The binning in sols also takes into
-! account the variability of Local Time in the bin and try to represent it too.
+! account the variability of Local Time in the bin and tries to represent it (with a number of 
+! LT values in each sol bin equal to the number of values in each MCS bin, and centered
+! around the corresponding MCS bin LT average).
 !
-! Eventually, the program outputs a netcdf file, filled with the output bin's data and edited with
+! Eventually, the program outputs a netcdf file, filled with the GCM binned data and edited with
 ! the same format than the MCS file (and with the same dimensions' declaration order).
 !
 ! Minimal requirements :
-!     - the MCS file must contain the variables dtemp,dtimeave,dtimemax,dtimemin,
-!                                               ntemp,ntimeave,ntimemax,ntimemin
-!     - the GCM file must : contain the variable temp (atmospheric temperature) ;
-!                           have the same altitude type as the MCS file ;
-!                           cover completely the NON-BINNED observation period
-!     - if the GCM file is a stats file, please name it with "stats" in it (to be changed in the future)
+!     - the MCS file must contain the variables dtemp,dtimeave,dtimemax,dtimemin,dnumbintemp,
+!                                               ntemp,ntimeave,ntimemax,ntimemin,nnumbintemp
+!     - the GCM file must : # contain the variable temp (atmospheric temperature) ;
+!                           # have the same altitude type as the MCS file ;
+!                           # cover completely the observation period
 !
 ! See also NOTA BENE in section 1.3
 ! More details can also be found in my internship report (in French) :
 !           /planeto/abierjon/rapportstage2019/Rapport_de_PFE_Master_Bierjon_2019.pdf
-! or directly in my office !
 !===================================================================================================
 
 
 use netcdf
+
 
 !===================================================================================================
@@ -43 +44 @@
 implicit none ! for no implicitly typed variables
 
+!------------------------
 ! Files:
 character(len=256) :: gcmfile ! GCM simulation file
@@ -49 +51 @@
 character(len=256) :: outfile ! output file
 
+!------------------------
 ! NetCDF stuff
 integer :: status ! NetCDF routines return code
@@ -55 +58 @@
 integer :: obsfid ! NetCDF observation file ID
 integer :: outfid ! NetCDF output file ID
-integer :: GCMvarid, OBSvarid, LT_id, outvarid ! to store the ID of a variable
+integer :: GCMvarid, OBSvarid, LT_id, numbin_id, outvarid ! to store the ID of a variable
 integer :: lat_dimid_obs,lon_dimid_obs,alt_dimid_obs,time_dimid_obs ! dimensions' ID in OBS and output files
 integer :: lat_dimid_gcm,lon_dimid_gcm,alt_dimid_gcm,time_dimid_gcm ! dimensions' ID in GCM file
-integer :: GCMvarshape(4), OBSvarshape(4), LTshape(3) ! to store a variable's coordinates order
-integer :: nb
-
+integer :: GCMvarshape(4), OBSvarshape(4), LTshape(3),numbinshape(4) ! to store a variable's coordinates order
+integer :: nb ! nb of dimensions of a variable
+
+!------------------------
+! Dimensions
 real,dimension(:),allocatable :: GCMlon, OBSlon ! longitude in the GCM & the Observer files
 integer GCMlonlen, OBSlonlen ! # of grid points along GCMlon & OBSlon
-
 real,dimension(:),allocatable :: GCMlat, OBSlat ! latitude in the GCM & the Observer files
 integer GCMlatlen, OBSlatlen ! # of grid points along GCMlat & OBSlat
-
 real,dimension(:),allocatable :: GCMalt, OBSalt ! can be geometric heights or pressure
 integer GCMaltlen, OBSaltlen ! # of grid point along GCMalt & OBSalt
 character (len=256) :: units ! to store the units attribute of altitude
 character(len=1) :: GCMalttype, OBSalttype ! altitude coord. type:'z' (altitude, m) 'p' (pressure, Pa)
-
 real,dimension(:),allocatable :: GCMtime, OBSLs ! time in the GCM diagfi (sols) & the Observer files (Ls)
 real,dimension(:),allocatable :: GCMstatstime ! time in the GCM stats file (LT at lon 0°)
@@ -77 +79 @@
 real :: starttimeoffset=0 ! offset (in sols) wrt Ls=0 of sol 0 in GCM file
 
+!------------------------
+! Variables
 character (len=64) :: GCMvarname, OBSvarname ! to store the name of the variable to retrieve
 real,dimension(:,:,:,:),allocatable :: GCM_var,OBS_var ! the 4D variable extracted from GCM & OBS files
 character (len=64) :: OBSLTave_name, OBSLTmax_name, OBSLTmin_name ! to store the name of the average, max and min of LT (ex : dtimeave, ntimemax...)
 real,dimension(:,:,:),allocatable :: OBSLT, OBSLTmax, OBSLTmin ! 3D variables extracted from obsfile (average, max and min of LT in a bin)
+character (len=64) :: numbin_name ! to store the name of the nb of values in an OBS temp bin (dnumbintemp/nnumbintemp)
+real,dimension(:,:,:,:),allocatable :: numbintemp ! nb of values in an OBS temp bin
 real :: GCMmiss_val, OBSmiss_val, LTmiss_val ! value to denote non-existant data
 real :: extr_value ! to store the result of the extraction subroutine
-
+real, dimension(:,:,:,:), allocatable :: outvar ! outvar(,,,): 4D array to store the output variable's data
+
+!------------------------
+! Time binning management
 real :: OBSdeltaLs ! difference of Ls between each observation bin
 real :: sol, maxsol, minsol ! sol date corresponding to Observed Ls and LT
 integer :: m_maxsol, m_minsol ! indexes of the maximum and minimum GCM sol taken in a bin for interpolation
 
-real,dimension(:,:,:),allocatable :: LTmod_arraymax, LTmod_arraymin ! declaration of the type of the function LTmod_array
+external LTmod ! declaration of the function LTmod
 real :: LTmod ! declaration of the type of the function LTmod
-real :: maxLTvar ! maximum variability of LT within one MCS bin
-integer :: LTcount ! nb of LT samples on which the interpolation is performed, for every LT interval
+integer :: LTcount ! nb of LT samples on which the interpolation is performed, for every LT interval = numbintemp[lon,lat,alt,Ls]
 
 integer :: solcount ! number of GCM sol integer values in one Ls interval
@@ -99 +107 @@
 integer :: errcount = 0 ! total number of GCM missing values
 real :: solbinned_value ! to store the extracted value averaged on sol samples, which is finally put in the output bin
-real :: GCMdeltat ! timestep of the GCM simulation (in Martian hours)
-
+
+!------------------------
+! Extraction & loop indexes
 real :: lon_val, lat_val, alt_val, Ls_val, LT_val ! where and when the output file is written at
 
@@ -106 +115 @@
 integer :: m,m_int,n ! sol binning loops iteration index
 
-real, dimension(:,:,:,:), allocatable :: outvar ! outvar(,,,): 4D array to store the output variable's data
 
 !===================================================================================================
@@ -228 +236 @@
 !===============================================================================
 
-! Observer variables to extract :
 !******************** NOTA BENE (cf sections 1.3 and 4)*************************
 ! We execute the program a first time with the daytime values, and then a second
@@ -234 +241 @@
 ! independent of this distinction, hence we execute them only in the first loop,
 ! when OBSvarname="dtemp". It also implies a lot of IF (when it's only executed in 
-! the first loop) or CASE (when it's executed in both loops but depends on wether 
+! the first loop) or CASE (when it's executed in both loops but depends on whether 
 ! it is daytime or nighttime).
-! This should be changed one day with a cleaner version of the code...
+! This may be changed one day with a cleaner version of the code...
 !*******************************************************************************
+
+! Observer variables to extract :
 OBSvarname="dtemp" ! we begin with daytime temperature
 write(*,*)"" ; write(*,*) "Beginning the 1st loop, on daytime values"; write(*,*)""
@@ -304 +313 @@
     OBSLTmax_name = "dtimemax"
     OBSLTmin_name = "dtimemin"
+    numbin_name   = "dnumbintemp"
   CASE ("ntemp")
     OBSLTave_name = "ntimeave"
     OBSLTmax_name = "ntimemax"
     OBSLTmin_name = "ntimemin"
+    numbin_name   = "nnumbintemp"
   END SELECT
   
+  ! Get the average values of LT in the OBS bins
   status=nf90_inq_varid(obsfid,trim(OBSLTave_name),LT_id)
-  error_text="Failed to find Observer local time ("//trim(OBSLTave_name)//") dimension"
+  error_text="Failed to find Observer local time ("//trim(OBSLTave_name)//") ID in "//trim(obsfile)
   call status_check(status,error_text)
   ! Sanity checks on the variable
@@ -351 +363 @@
   ! Get the maximum values of LT in the OBS bins
   status=nf90_inq_varid(obsfid,trim(OBSLTmax_name),LT_id)
-  error_text="Failed to find Observer max local time ("//trim(OBSLTmax_name)//") dimension"
+  error_text="Failed to find Observer max local time ("//trim(OBSLTmax_name)//") ID in "//trim(obsfile)
   call status_check(status,error_text)
   ! Sanity checks on the variable
@@ -386 +398 @@
   ! Get the minimum values of LT in the OBS bins
   status=nf90_inq_varid(obsfid,trim(OBSLTmin_name),LT_id)
-  error_text="Failed to find Observer min local time ("//trim(OBSLTmin_name)//") dimension"
+  error_text="Failed to find Observer min local time ("//trim(OBSLTmin_name)//") ID in "//trim(obsfile)
   call status_check(status,error_text)
   ! Sanity checks on the variable
@@ -418 +430 @@
   write(*,*) "Loaded ", trim(OBSLTmin_name), " from the obsfile"
 
-  
-  SELECT CASE (OBSvarname)
-  CASE ("dtemp")
-    maxLTvar=maxval((OBSLTmax(:,:,:)-OBSLTmin(:,:,:)))
-    write(*,*)"The max daytime LT variability is ",maxLTvar,"hours (within a same bin)"
-  CASE ("ntemp")
-    ! in case of night local times, we replace hours from a [0;24[ interval to a 
-    ! [-12;12[ interval in which midnight = 0, with the subroutine LTmod
-    allocate(LTmod_arraymax(OBSlonlen,OBSlatlen,OBSLslen))
-    allocate(LTmod_arraymin(OBSlonlen,OBSlatlen,OBSLslen))
-    do i=1,OBSlonlen
-      do j=1,OBSlatlen
-        do l=1,OBSLslen
-          LTmod_arraymax(i,j,l) = LTmod(OBSLTmax(i,j,l))
-          LTmod_arraymin(i,j,l) = LTmod(OBSLTmin(i,j,l))
-        enddo
-      enddo
-    enddo
-    maxLTvar=maxval(LTmod_arraymax(:,:,:) - LTmod_arraymin(:,:,:))
-    write(*,*)"The max nighttime LT variability is ",maxLTvar,"hours (within a same bin)"
-  END SELECT
-  
-  LTcount=floor(maxLTvar) ! integer nb of hours fitting in the broadest interval of LT
+
+  ! Get the number of values of temp in the OBS bins
+  status=nf90_inq_varid(obsfid,trim(numbin_name),numbin_id)
+  error_text="Failed to find Observer nb of temp values in bin ("//trim(numbin_name)//")'s ID in "//trim(obsfile)
+  call status_check(status,error_text)
+  ! Sanity checks on the variable
+  status=nf90_inquire_variable(obsfid,numbin_id,ndims=nb,dimids=numbinshape)
+  error_text="Failed to obtain information on variable "//trim(numbin_name)
+  call status_check(status,error_text)
+  ! Check that it is a 4D variable
+  if (nb.ne.4) then
+    write(*,*) "Error, expected a 4D (lon-lat-alt-time) variable for of ", trim(numbin_name), "!"
+    stop
+  endif
+  ! Check that its dimensions are indeed lon,lat,alt,time (in the right order)
+  if (numbinshape(1).ne.lon_dimid_obs) then
+    write(*,*) "Error, expected first dimension of ", trim(numbin_name), " to be longitude!"
+    stop
+  endif
+  if (numbinshape(2).ne.lat_dimid_obs) then
+    write(*,*) "Error, expected second dimension of ", trim(numbin_name), " to be latitude!"
+    stop
+  endif
+  if (numbinshape(3).ne.alt_dimid_obs) then
+    write(*,*) "Error, expected third dimension of ", trim(numbin_name), " to be altitude!"
+    stop
+  endif
+  if (numbinshape(4).ne.time_dimid_obs) then
+    write(*,*) "Error, expected fourth dimension of ", trim(numbin_name), " to be time!"
+    stop
+  endif
+  ! Length allocation for each dimension of the 4D variable
+  allocate(numbintemp(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
+  ! Load datasets
+  status=nf90_get_var(obsfid,numbin_id,numbintemp)
+  error_text="Failed to load "//trim(numbin_name)//" from the obsfile"
+  call status_check(status,error_text)
+  write(*,*) "Loaded ", trim(numbin_name), " from the obsfile"
 
 !===================================================================================================
@@ -453 +480 @@
     write(*,*)"";write(*,*) "-> Enter input file name (GCM simulation) :"
     read(*,*) gcmfile
-    if (index(gcmfile,"stats").ne.0) then ! if the GCM file name contains "stats", we assume it is a stats file
-    ! => A REMPLACER PAR UN TEST SUR LA DIMENSION TIME ? REGARDER DANS AUTRES UTILITAIRES 
-      is_stats = .true.
-      write(*,*)"The GCM file is recognized as a stats file."
-    else
-      write(*,*)"The GCM file is recognized as a diagfi/concatnc file."
-    endif
 
     ! Open GCM file
@@ -512 +532 @@
     ! GCM Time
     status=nf90_inq_dimid(gcmfid,"Time",time_dimid_gcm)
-    error_text="Failed to find GCM time (sols) dimension"
-    call status_check(status,error_text)
-
-    if (.not.is_stats) then 
-      status=nf90_inquire_dimension(gcmfid,time_dimid_gcm,len=GCMtimelen)
-      error_text="Failed to find GCM time length"
-      call status_check(status,error_text)
-      allocate(GCMtime(GCMtimelen))
-    else 
-      status=nf90_inquire_dimension(gcmfid,time_dimid_gcm,len=GCMstatstimelen)
-      error_text="Failed to find GCM stats time length"
-      call status_check(status,error_text)
-      allocate(GCMstatstime(GCMstatstimelen))
-    endif
-
+    error_text="Failed to find GCM time dimension"
+    call status_check(status,error_text)
+
+    status=nf90_inquire_dimension(gcmfid,time_dimid_gcm,len=GCMtimelen)
+    error_text="Failed to find GCM time length"
+    call status_check(status,error_text)
+    allocate(GCMtime(GCMtimelen))
+    
     status=nf90_inq_varid(gcmfid,"Time",GCMvarid)
     error_text="Failed to find GCM time ID"
     call status_check(status,error_text)
-
-    ! Read GCMtime
-    if (.not.is_stats) then ! if GCM file is a diagfi, time is in sols at longitude 0°
-      status=nf90_get_var(gcmfid,GCMvarid,GCMtime)
-      error_text="Failed to load GCMtime (sols)"
-      call status_check(status,error_text)
-    else ! if GCM file is a stats, time is in LT at longitude 0°
+    
+    status=nf90_get_var(gcmfid,GCMvarid,GCMtime)
+    error_text="Failed to load GCMtime"
+    call status_check(status,error_text)
+
+    ! is_stats ?
+    if ((GCMtimelen.eq.12).and.(GCMtime(1).eq.2.).and.(GCMtime(GCMtimelen).eq.24.)) then 
+      ! if GCM file is a stats, time is in LT at longitude 0° and not in sols at longitude 0°
+      write(*,*)"The GCM file is recognized as a stats file."
+      is_stats = .true.
+      deallocate(GCMtime)
+      GCMstatstimelen = GCMtimelen
+      allocate(GCMstatstime(GCMstatstimelen))
       status=nf90_get_var(gcmfid,GCMvarid,GCMstatstime)
       error_text="Failed to load GCMstatstime (LT at lon 0)"
       call status_check(status,error_text)
+    else
+      write(*,*)"The GCM file is recognized as a diagfi/concatnc file."
     endif
 
@@ -550 +571 @@
       GCMtime(:)=GCMtime(:)+starttimeoffset
     endif 
-
-    ! Get the timestep of the simulation (in Martian hours)
-    if (.not.is_stats) then
-      GCMdeltat = (GCMtime(2)-GCMtime(1))*24.
-    else
-      GCMdeltat = (GCMstatstime(2)-GCMstatstime(1))
-    endif
-    write(*,*)"GCM simulation's timestep is ", GCMdeltat, " hours."
 
     ! Check of temporal coherence between gcmfile & obsfile
@@ -749 +762 @@
   error_text="Error: could not define the variable "//trim(OBSLTave_name)//" in the outfile"
   call status_check(status,error_text)
+  
+  status=nf90_def_var(outfid,trim(numbin_name),nf90_float,numbinshape,numbin_id)
+  error_text="Error: could not define the variable "//trim(numbin_name)//" in the outfile"
+  call status_check(status,error_text)
 
   ! Write the attributes
@@ -792 +809 @@
   write(*,*)"Local Time (", trim(OBSLTave_name), ") has been created in the outfile"
   write(*,'("  with missing_value attribute : ",1pe12.5)')OBSmiss_val
+  
+  !! Number of temp values in bin
+  SELECT CASE (OBSvarname)
+  CASE ("dtemp")
+    status=nf90_put_att(outfid,numbin_id,"long_name","Number of temp values in bin - day side")
+  CASE ("ntemp")
+    status=nf90_put_att(outfid,numbin_id,"long_name","Number of temp values in bin - night side")
+  END SELECT
+  error_text="Error: could not put the long_name attribute for the variable "//trim(numbin_name)//" in the outfile"
+  call status_check(status,error_text)
+
+  write(*,*)"Number of temp values in bin (", trim(numbin_name), ") has been created in the outfile"
+  
 
   ! End the netcdf define mode (and thus enter the "data writing" mode)
@@ -867 +897 @@
         if (lon_val.gt.180.) lon_val=lon_val-360.
 
-        LT_val=OBSLT(i,j,l) ! find the Observer corresponding LT value at lon_val, lat_val, Ls_val
-        
+        LT_val=OBSLT(i,j,l) ! find the Observer average LT value at bin(lon_val, lat_val, Ls_val)
+        
         if ((LT_val.lt.0.).or.(LT_val.gt.24.)) then
           if (LT_val.eq.LTmiss_val) then
@@ -880 +910 @@
           endif
         endif
-        ! Generate the sol list for the interpolation
-        allocate(sol_list(solcount*LTcount))
-        call gen_sol_list(solcount,int_sol_list,LTcount,LT_val,OBSLTmax(i,j,l),OBSLTmin(i,j,l),lon_val,LTmod,OBSvarname,&
-                          sol_list)
 
         altitude: do k=1,OBSaltlen ! loop on all observed altitudes
@@ -897 +923 @@
 !! 3.3.2 Compute GCM sol date corresponding to Observer Ls 
 !!       (via m_(min/max)sol) and LT (via OBSLT(min/max))
-!!===============================================================       
+!!===============================================================
+          LTcount=floor(numbintemp(i,j,k,l)) ! find the Observer number of temp values at bin(lon_val, lat_val, alt_val, Ls_val)
+          if (LTcount.eq.0.) then
+            ! Write a missing_value to output
+            outvar(i,j,k,l)=OBSmiss_val
+            CYCLE altitude ! go directly to the next altitude
+          endif
+          if (LTcount.lt.0.) then
+            write(*,*) "Unexpected value for LTcount: ",LTcount
+            stop
+          endif
+
+          ! Generate the sol list for the interpolation
+          allocate(sol_list(solcount*LTcount))
+          call gen_sol_list(solcount,int_sol_list,LTcount,LT_val,OBSLTmax(i,j,l),OBSLTmin(i,j,l),lon_val,LTmod,OBSvarname,&
+                            sol_list) 
+                
           solerrcount=0
           solbinned_value=0
@@ -929 +971 @@
 
           errcount=errcount+solerrcount
-
+          
+          deallocate(sol_list)
+          
         enddo altitude ! end loop on observed altitudes
-
-        deallocate(sol_list)
-
       enddo longitude ! end loop on observed longitudes
     enddo latitude ! end loop on observed latitudes
@@ -951 +992 @@
 
   status = nf90_put_var(outfid, LT_id, OBSLT) ! write the MCS d/ntimeave as the output Local Time
-  error_text="Error: could not write Local Time data in the outfile"
-  call status_check(status,error_text)
-
+  error_text="Error: could not write "//trim(OBSLTave_name)//" data in the outfile"
+  call status_check(status,error_text)
+  
+  status = nf90_put_var(outfid, numbin_id, numbintemp)
+  error_text="Error: could not write "//trim(numbin_name)//" data in the outfile"
+  call status_check(status,error_text)
+  
 !!==============================================================================  
 !! 4. End of the Day/Night loop
@@ -965 +1010 @@
     deallocate(OBSLTmax)
     deallocate(OBSLTmin)
+    deallocate(numbintemp)
     deallocate(outvar)
     write(*,*)"";write(*,*)""; write(*,*) "Beginning the 2nd loop, on nighttime values" ; write(*,*)""
@@ -1562 +1608 @@
 real, intent(in) :: lon ! longitude value for the transformation into a sol value at lon=0°
 external f_LT ! function that changes LT interval for night LT
-real :: f_LT
+real f_LT
 character (len=64), intent(in) :: LTname ! to know if we have day or night values
 
@@ -1657 +1703 @@
 ! night local times)
 !================================================================
-real :: LT
+implicit none
+real,intent(in) :: LT
 
 LTmod = 2*mod(LT,12.)-LT