Changeset 2308 for trunk/LMDZ.MARS/util/concatnc.F90

Timestamp:

May 4, 2020, 11:07:20 PM (5 years ago)

Author:

abierjon

Message:

Mars GCM:
Following r2303 (concatenate concat files) and truly resolving Ticket #46 :
1) little bugs fixed from the previous revision ; 2) add the possibility to change the
offset of the output file time axis by asking the user the new offset and the level of
priority to put on it over the starting day stored in the first input file
(NB: this last point also implies a change in concatnc.def, but retrocompatibility
with old concatnc.def files has been ensured)

AB

File:

• trunk/LMDZ.MARS/util/concatnc.F90 (modified) (14 diffs)

trunk/LMDZ.MARS/util/concatnc.F90

@@ -20 +20 @@
 ! + read and write controle field, if available. TN, October 2013
 ! + possibility to concatenate concat files with a coherent 
-!   time axis in the output. AB, April 2020
+!   time axis in the output, and to redistribute the time axis
+!   from an offset given by the user. AB, April 2020
 ! ********************************************************
 
@@ -116 +117 @@
 real :: memotime
 ! memotime: (cumulative) time value, in martian days (sols)
+real :: starttimeoffset
+! starttimeoffset: offset given by the user for the output time axis
+character (len=1) :: firstsol
+! firstsol: =[y/n], to know whether the user offset or the first file's starting day must be used in the output
 real :: missing
 !PARAMETER(missing=1E+20)
@@ -126 +131 @@
 !==============================================================================
 write(*,*) 
-write(*,*) "-- Program written specifically for NetCDF 'diagfi'/'concat' files from the martian GCM --"
+write(*,*) "-- Program written specifically for NetCDF 'diagfi'/'concat' files"
+write(*,*) "                     from the martian GCM --"
 write(*,*) 
 write(*,*) "which files do you want to use?"
@@ -145 +151 @@
 
 !==============================================================================
-! 1.2. Ask for starting day value (memotime)
+! 1.2. Ask for starting day value (starttimeoffset)
 !==============================================================================
 
 write(*,*)
-!write(*,*) "Beginning day of the first specified file?"
-write(*,*) "Starting day of the run stored in the first input file?"
-write(*,*) " (Obsolete if the controle field is present, answer any number)"
-write(*,*) "(e.g.: 100 if that run started at time=100 sols)"
-read(*,*) memotime
+write(*,*) "Starting day to be stored in the output file time axis?"
+write(*,*) "(e.g.: 100 if you want the output file to start at time=100 sols)"
+read(*,*) starttimeoffset
 
 !==============================================================================
@@ -174 +178 @@
 !==============================================================================
 
+write(*,*)
 write(*,*) "Warning: to read the result with grads, choose sol"
 write(*,*) "Warning: use ferret to read the non linear scale ls"
@@ -322 +327 @@
 ! 2.2. Read (and check) dimensions of variables from input file
 !==============================================================================
-
+!==============================================================================
+! 2.2.1 Lat, lon, alt, GCM_layers
+!==============================================================================
    ierr=NF_INQ_DIMID(nid,"latitude",latdim)
    ierr=NF_INQ_VARID(nid,"latitude",latvar)
@@ -349 +356 @@
    ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
 !  write(*,*) "altlen: ",altlen
-
-   ierr=NF_INQ_DIMID(nid,"index",ctldim)
-   ierr=NF_INQ_VARID(nid,"controle",ctlvar)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'Field <controle> is missing in file '//file(i)
-      ctllen=0
-      !stop ""
-   else
-      ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
-   endif
-!  write(*,*) "controle: ",controle
 
 ! load size of aps() or sigma() (in case it is not altlen)
@@ -375 +371 @@
 
 !==============================================================================
+! 2.2.1 Check presence and read "controle" dimension & variable from input file
+!==============================================================================
+   ierr=NF_INQ_DIMID(nid,"index",ctldim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'Dimension <index> is missing in file '//file(i)
+      write(*,*) "The program continues..."
+      ctllen=0
+      !stop ""
+   else
+      ierr=NF_INQ_VARID(nid,"controle",ctlvar)
+      if (ierr.NE.NF_NOERR) then
+         write(*,*) 'Field <controle> is missing in file '//file(i)
+         write(*,*) "The program continues..."
+         ctllen=0
+         !stop ""
+      else
+         ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
+      endif
+   endif
+!  write(*,*) "controle: ",controle
+
+   if (ctllen .ne. 0) then ! variable controle
+      allocate(ctl(ctllen))
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
+#else
+      ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
+#endif
+      ctlsol = ctl(4)
+   endif
+
+!==============================================================================
+! 2.3 Read "Time" dimension & variable from input file
+!==============================================================================
+   ierr=NF_INQ_DIMID(nid,"Time",timedim)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Dimension <Time> is missing in file '//file(i)
+      stop ""
+   endif
+   ierr=NF_INQ_VARID(nid,"Time",timevar)
+   if (ierr.NE.NF_NOERR) then
+      write(*,*) 'ERROR: Field <Time> is missing in file '//file(i)
+      stop ""
+   endif
+   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
+!  write(*,*) "timelen: ",timelen
+
+   ! allocate time() array and fill it with values from input file
+   allocate(time(timelen))
+
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
+#else
+   ierr = NF_GET_VAR_REAL(nid,timevar,time)
+#endif
+
+
+!==============================================================================
 ! 2.3. Read (and check compatibility of) dimensions of
 !       variables from input file
 !==============================================================================
-
-   if (ctllen .ne. 0) then ! variable controle
-      allocate(ctl(ctllen))
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
-#else
-      ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
-#endif
-      ctlsol = ctl(4)
-   endif
 
    if (i==1) then ! First call; initialize/allocate
@@ -406 +450 @@
       ierr = NF_GET_VAR_REAL(nid,altvar,alt)
 #endif
+      firstsol = 'n' ! defaut value
+      write(*,*)
       if (ctllen .ne. 0) then
-         if (modulo(int(memotime),669)/=modulo(int(ctlsol),669)) then
-           write(*,*) "WARNING: Starting day of the run is not ",&
-                                modulo(int(memotime),669)," but ",modulo(int(ctlsol),669),"!!"
-           write(*,*) "Starting day of the run has been corrected."
-!           memotime=float(modulo(int(ctl(4)),669)) + ctl(27)
-           memotime=float(int(ctlsol)) + ctl(27)
-           ctl(4) = 0.  ! values written in the output
-           ctl(27) = 0. ! file controle variable (resp. day_ini and hour_ini)
+         if (int(starttimeoffset) .ne. (int(ctlsol)+int(time(1))) ) then
+           write(*,*) "WARNING: Starting day of the first file is not ",&
+                                int(starttimeoffset)," but ",int(ctlsol)+int(time(1)),"!"
+           write(*,*) "Answer:"
+           write(*,*) " y to use the starting day you have inputted"
+           write(*,*) " n to keep the starting day of the first file"
+           write(*,*) "as the starting day of the output file"
+           read(*,*) firstsol
          endif
+         
+         if (firstsol.eq.'y') then
+           write(*,*) "The day given by the user is used as the starting day of the ouput file"
+           starttimeoffset = float(int(ctlsol)+int(time(1))) + ctl(27) - starttimeoffset
+         else ! if firstsol.eq.'n' or other (default case, for retrocompatibility with previous concatnc.def)
+           write(*,*) "The starting day of the 1st file is used as the starting day of the ouput file"
+           starttimeoffset = 0
+         endif
+         
+         memotime=float(int(ctlsol)+int(time(1))) + ctl(27)
+         ctl(4) = 0.  ! values written in the output
+         ctl(27) = 0. ! file controle variable (resp. day_ini and hour_ini)
+         
+      else ! if no variable "controle"
+         if (int(starttimeoffset) .ne. (int(time(1))) ) then
+           write(*,*) "WARNING: Starting day of the first file is not ",&
+                                int(starttimeoffset)," but ",int(time(1)),"!"
+           write(*,*) "Answer:"
+           write(*,*) " y to use the starting day you have inputted"
+           write(*,*) " n to keep the starting day of the first file"
+           write(*,*) "as the starting day of the output file"
+           read(*,*) firstsol
+         endif
+         
+         if (firstsol.eq.'y') then
+           write(*,*) "The day given by the user is used as the starting day of the ouput file"
+           starttimeoffset = int(time(1)) - starttimeoffset
+         else ! if firstsol.eq.'n' or other (default case, for retrocompatibility with previous concatnc.def)
+           write(*,*) "The starting day of the 1st file is used as the starting day of the ouput file"
+           starttimeoffset = 0
+         endif
+
+         memotime=float(int(time(1)))
       endif
+        
    ! Initialize output file's lat,lon,alt and time dimensions
+      write(*,*)
       call initiate (filename,lat,lon,alt,ctl,GCM_layers,nout,&
        latdimout,londimout,altdimout,timedimout,&
@@ -429 +510 @@
    ! Check that latitude,longitude and altitude of current input file
    ! are identical to those of the output file
+   ! and that either all the files or no file contain the controle variable
       if (memolatlen/=latlen) then
            write(*,*) "ERROR: Not the same latitude axis"
@@ -444 +526 @@
    endif ! of if (i==1)
 
-!==============================================================================
-! 2.4. Handle "Time" dimension from input file
-!==============================================================================
-
-!==============================================================================
-! 2.4.1 Read "Time" dimension from input file
-!==============================================================================
-   ierr=NF_INQ_DIMID(nid,"Time",timedim)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Dimension <Time> is missing in file '//file(i)
-      stop ""
-   endif
-   ierr=NF_INQ_VARID(nid,"Time",timevar)
-   if (ierr.NE.NF_NOERR) then
-      write(*,*) 'ERROR: Field <Time> is missing in file '//file(i)
-      stop ""
-   endif
-   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
-!  write(*,*) "timelen: ",timelen
-
-   ! allocate time() array and fill it with values from input file
-   allocate(time(timelen))
-
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
-#else
-   ierr = NF_GET_VAR_REAL(nid,timevar,time)
-#endif
-
-!==============================================================================
-! 2.4.2 Write/extend "Time" dimension/values in output file
+
+!==============================================================================
+! 2.4 Write/extend "Time" dimension/values in output file
 !==============================================================================
 
    rep=0
    write(*,*)
-   write(*,'(a3,1x,f6.1)') "Sol",memotime
+   write(*,'(a3,1x,f6.1)') "Sol",memotime-starttimeoffset
 !   if (ctllen.ne.0) write(*,'(a30,1x,f6.1)') "In the current file's ctl: Sol",ctlsol
 
@@ -486 +540 @@
      do k=reptime+1,reptime+timelen
        rep=rep+1
-       if ((time(rep)+ctlsol).le.memotime)then
+       if ((time(rep)+ctlsol).le.memotime) then
          write(*,*) "ERROR : the time intervals of the files are not dissociated"
          stop ""
        endif
 #ifdef NC_DOUBLE
-       ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,k,1,ctlsol+time(rep))
-#else
-       ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,ctlsol+time(rep))
-#endif
+       ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,k,1,ctlsol+time(rep)-starttimeoffset)
+#else
+       ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,ctlsol+time(rep)-starttimeoffset)
+#endif
+       ! 
      enddo
 
-   else
+   ! Compute new time offset (for further concatenations)
+   memotime=ctlsol+time(timelen) ! the actual last day of the current file
+                                 ! becomes the reference time for next file
+                                 ! leaving potential "time holes" in the output time axis
+                                 ! at the junction of two files' time axes
+   
+   else ! if no variable "controle"
+   
      do k=reptime+1,reptime+timelen
         rep=rep+1
 #ifdef NC_DOUBLE
-        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,k,1,memotime+time(rep))
-#else
-        ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,memotime+time(rep))
+        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,k,1,memotime+time(rep)-int(time(1))-starttimeoffset)
+#else
+        ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,memotime+time(rep)-int(time(1))-starttimeoffset)
 #endif
      enddo
+     
+   ! Compute new time offset (for further concatenations)
+   memotime=memotime+time(timelen)-time(1) ! files are concatenated one after another
+                                           ! even if their time axes normally don't 
+                                           ! directly follow each other
    endif ! if ctllen.ne.0
-
-   ! Compute new time offset (for further concatenations)
-   memotime=memotime+time(timelen)
 
 !==============================================================================
@@ -648 +712 @@
 
 if (axis=="ls") then
+   write(*,*);write(*,*) "Converting the time axis from sol to Ls..."
    call change_time_axis(nout,ierr)
 endif
@@ -654 +719 @@
 ierr=NF_CLOSE(nout)
 
+write(*,*);write(*,*) "Program concatnc completed !"
 contains