Changeset 2450 for trunk/LMDZ.MARS/util

Timestamp:

Jan 25, 2021, 8:03:25 PM (4 years ago)

Author:

abierjon

Message:

Mars GCM:
Changes in lslin.F90 :

• make sure the interpolation takes into account potential missing values in the input file fields, especially when the average/binning option is not used (ticket #68)
• enable the averaging option even when the Ls timestep is set automatically /!\ previous lslin.def must be updated in consequence !
• code clean-up, especially remove duplicates on controle field and altitude units
• replace all the tests "if (ndim.ge.3) then" by a unique test at the beginning of the var loop
• change "title" attribute into "long_name" by default (ticket #48)

+ add an indicative updated lslin.def
+ update util/README

AB

Location:

trunk/LMDZ.MARS/util

Files:

• README (modified) (2 diffs)
• lslin.F90 (modified) (28 diffs)
• lslin.def (added)

trunk/LMDZ.MARS/util/README

@@ -90 +90 @@
 
 --------------------------------------------------------------------
-4) lslin.e : redistribute and AVERAGE gcm output in Ls coordinate.
---------------------------------------------------------------------
-
-This program has been designed to interpol data in Solar Longitude (Ls)
+4) lslin.e : redistribute and average gcm output in Ls coordinate.
+--------------------------------------------------------------------
+
+This program has been designed to interpolate data in Solar Longitude (Ls)
 linear time coordinate (usable with grads) from Netcdf diagfi or concatnc  
-files and to AVERAGE the instantaneous data in Ls period (for comparison with
-binned dataset, for instance).
+files, and to average the instantaneous data in periodic Ls intervals
+(for comparison with binned dataset, for instance).
 
 lslin also creates a lslin.ctl file that can be read
@@ -103 +103 @@
 is too small"...
 
-Output file is : name_of_input_file_LS.nc
+Output file is : name_of_input_file_Ls.nc
 
 --------------------------------------------------------------------

trunk/LMDZ.MARS/util/lslin.F90

@@ -23 +23 @@
 character (len=50), dimension(:), allocatable :: var
 ! var(): names of the variables
-character (len=50) :: title,units
-! title(),units(): [netcdf] "title" and "units" attributes
-character (len=50) :: units_alt
-! var(): units of altitude which can be 'km' or 'Pa' (after zrecast)
+character (len=50) :: long_name,units
+! long_name(),units(): [netcdf] "long_name" and "units" attributes
 character (len=50) :: altlong_name,altunits,altpositive
 ! altlong_name(): [netcdf] altitude "long_name" attribute
@@ -37 +35 @@
 ! vartmp(): used for temporary storage of various strings
 character (len=1) :: answer
-! answer: the character 'y' (or 'Y'), if input file is of 'concatnc' type
+! answer: user's answers at multiple questions (especially, answer='y'/'Y')
 character (len=1) :: average
 ! average: the character 'y' to average within the Ls timestep
@@ -55 +53 @@
 integer :: day_ini ! Ehouarn: integer or real ?
 ! day_ini: starting date/time of the run stored in input file
-integer, parameter :: length=100
-! length: (default) lenght of tab_cntrl()
-real, dimension(length) :: tab_cntrl
-! tab_cntrl(length): array, stored in the input file,
-!                which contains various control parameters.
 real, dimension(:), allocatable:: lat,lon,alt,time,lsgcm,timels,ctl
 ! lat(): latitude coordinates (read from input file)
@@ -129 +122 @@
 read(*,'(a50)') infile 
 
-!write(*,*) "Is it a concatnc file? (y/n)?"
-write(*,*) "Do you want to specify the beginning day of the file"
-write(*,*) "in case the controle field is not present ? (y/n)?"
-read(*,*) answer
-if ((answer=="y").or.(answer=="Y")) then
-   write(*,*) "Beginning day of the file?"
-   read(*,*) reptime 
-!   start_var=8 ! 'concatnc' type of file
-!else
-!   start_var=16 ! 'diagfi' type of file
-! N.B.: start_var is now computed; see below
-endif
-
-
 !==============================================================================
 ! 1.2. Open input file and read/list the variables it contains
@@ -178 +157 @@
 ierr=NF_INQ_NVARS(nid,nbvarfile)
 ! nbvarfile is now equal to the (total) number of variables in input file
-
-allocate(var(nbvarfile-(start_var-1)))
-
-! Yield the list of variables (to the screen)
-write(*,*)
-do i=start_var,nbvarfile
-   ierr=NF_INQ_VARNAME(nid,i,vartmp)
-   write(*,*) trim(vartmp)
-enddo
-!nbvar=0
-
-! Ehouarn: Redundant (and obsolete) lines ?
-!  nbvar=nbvarfile-(start_var-1)
-!  do j=start_var,nbvarfile
-!  ierr=nf_inq_varname(nid,j,var(j-(start_var-1)))
-!  enddo
                                        
 ! Get the variables' names from the input file (and store them in var())
+write(*,*) ""
 nbvar=nbvarfile-(start_var-1)
+allocate(var(nbvar))
 do i=1,nbvar
    ierr=NF_INQ_VARNAME(nid,i+(start_var-1),var(i))
@@ -205 +170 @@
 ! 1.3. Output file name
 !==============================================================================
-! outfile="lslin.nc"
 outfile=infile(1:len_trim(infile)-3)//"_Ls.nc"
- write(*,*) outfile
+write(*,*) ""
+write(*,*) "Output file : "//trim(outfile)
 
 
@@ -247 +212 @@
    else
       ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
-      units_alt=""
-      ierr=NF_GET_ATT_TEXT(nid,varid,"units",units_alt)
-   endif
-   write(*,*) "altlen: ",altlen
-! Get altitude attributes to handle files with any altitude type
-   ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
-   ierr=nf_get_att_text(nid,altvar,'units',altunits)
-   ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
+      
+    ! Get altitude attributes to handle files with any altitude type
+      ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
+      ierr=nf_get_att_text(nid,altvar,'units',altunits)
+      ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
+   endif
+!   write(*,*) "altlen: ",altlen
 
 
 ! Allocate lat(), lon() and alt()
-      allocate(lat(latlen))
-      allocate(lon(lonlen))
-      allocate(alt(altlen))
+   allocate(lat(latlen))
+   allocate(lon(lonlen))
+   allocate(alt(altlen))
 
 ! Read lat(),lon() and alt() from input file
-      ierr = NF_GET_VAR_REAL(nid,latvar,lat)
-      ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
-      ierr = NF_GET_VAR_REAL(nid,altvar,alt)
-      if (ierr.NE.NF_NOERR) then
-          if (altlen.eq.1) alt(1)=0.
-      end if
+   ierr = NF_GET_VAR_REAL(nid,latvar,lat)
+   ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
+   ierr = NF_GET_VAR_REAL(nid,altvar,alt)
+   if (ierr.NE.NF_NOERR) then
+       if (altlen.eq.1) alt(1)=0.
+   end if
 
 
@@ -304 +268 @@
 
 !==============================================================================
-! 2.2. Create (and initialize latitude, longitude and altitude in) output file
+! 2.2. Create (and initialize) latitude, longitude and altitude in output file
 !==============================================================================
 
@@ -316 +280 @@
               nout,londimout,latdimout,altdimout)
 
+   write(*,*)""
 !==============================================================================
 ! 2.3. Read time from input file
@@ -339 +304 @@
 allocate(lsgcm(timelen))
 
-
-
-
 ierr = NF_GET_VAR_REAL(nid,timevar,time)
 
-
 !==============================================================================
 ! 2.4. Initialize day_ini (starting day of the run)
 !==============================================================================
 
-
-write(*,*) 'answer',answer
+if (ctllen .ne. 0) then ! controle is present in input file
+  write(*,*) "Do you want to specify the beginning day of the file (y/n) ?"
+  write(*,*) "If 'n', the value stored in the controle field will be used."
+  read(*,*) answer
+else ! controle is not present in the input file
+  write(*,*) "There is no field controle in file "//trim(infile)
+  write(*,*) "You have to specify the beginning day of the file."
+  answer="y"
+endif
+
+if ((answer=="y").or.(answer=="Y")) then
+   write(*,*) "Beginning day of the file?"
+   read(*,*) reptime
+endif
+
+write(*,*) 'answer : ',answer
 if ((answer/="y").and.(answer/="Y")) then
-   ! input file is of 'concatnc' type; the starting date of the run
+   ! the starting date of the run
    ! is stored in the "controle" array
-   ierr = NF_INQ_VARID (nid, "controle", varid)
-   if (ierr .NE. NF_NOERR) then
-      write(*,*) 'ERROR: <controle> variable is missing'
-      write(*,*) NF_STRERROR(ierr)
-      stop ""
-   endif
-
-
-
-
-   ierr = NF_GET_VAR_REAL(nid, varid, tab_cntrl)
-
-
-   if (ierr .NE. NF_NOERR) then
-      write(*,*) 'ERROR: Failed to read <controle>'
-      write(*,*) NF_STRERROR(ierr)
-     stop ""
-   endif
-
-   day_ini = nint(tab_cntrl(4))                                               
+   
+   day_ini = nint(ctl(4))                                               
    day_ini = modulo(day_ini,669)                                                    
    write(*,*) 'day_ini', day_ini
@@ -397 +354 @@
 write(*,*) 'Input data LS range :', lsgcm(1),lsgcm(timelen)
 
-!******************************************
-write(*,*) "Automatic Ls timetsep (y/n)?"
+! generate the time step of the new Ls axis for the output
+write(*,*) "Automatic Ls timestep (y/n)?"
 read(*,*) answer
 if ((answer=="y").or.(answer=="Y")) then
@@ -417 +374 @@
     ls0=lsgcm(1) ! First Ls date
 else
-    write(*,*) "New timestep in Ls (deg)"
+    write(*,*) "New timestep in Ls (deg) :"
     read(*,*) deltals
-    write(*,*) "First Ls (deg)"
+    write(*,*) "First Ls (deg) :"
     read(*,*) ls0
     if (ls0.lt.lsgcm(1)) then 
@@ -426 +383 @@
     end if
     write(*,*) "First Ls date (deg) = ", ls0  ! FF2013
-
-    do while((average.ne.'y').and.(average.ne.'n')) 
-       write(*,*) "Average data within the Ls timestep (y/n) or interpolate ? "
-       read(*,*) average
-    end do
-endif
-  NLs=int(Int((lsgcm(timelen)-ls0)/deltals) +1)   ! FF2013
-!********************************************
+endif
+
+NLs=int(int((lsgcm(timelen)-ls0)/deltals) +1)   ! FF2013
 allocate(timels(Nls))
 
 ! Build timels()
-timels(1) = 0.01*nint(100.*ls0) 
+timels(1) = 0.01*nint(100.*ls0) ! 2 decimals
 do k=2,Nls
    timels(k) = timels(k-1) + deltals
 end do
 
-write(*,*) 'timestep in Ls (deg) ', deltals
-write(*,*) 'output data LS range : ', timels(1),timels(Nls)
+write(*,*) 'New timestep in Ls (deg) ', deltals
+write(*,*) 'Output data LS range : ', timels(1),timels(Nls)
+
+! create averaged bins or interpolate at exact Ls ?
+write(*,*) "Average data within the Ls timestep (y/n) or interpolate ? "
+read(*,*) average
+write(*,*)""
 
 !==============================================================================
@@ -450 +407 @@
 
 do k=1,Nls
-
-
-
-ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,timels(k))
-
+  ierr= NF_PUT_VARA_REAL(nout,timevarout,k,1,timels(k))
 enddo
-
 
 !==============================================================================
@@ -474 +426 @@
       write(*,*) 'ERROR: Field <',var(j),'> not found'
       write(*,*) NF_STRERROR(ierr)
-      stop ""
-   endif
+      stop "Better stop now..."
+   endif
+   
    ! Get the value of 'ndim' for this varriable
    ierr=NF_INQ_VARNDIMS(nid,varid,ndim)
    write(*,*) 'ndim', ndim
+   
+   ! Check that it is a 3D or 4D variable
+   if (ndim.lt.3) then
+     write(*,*) "Error:",trim(var(j))," is neither a 3D nor a 4D variable"
+     write(*,*) "We'll skip that variable..."
+     CYCLE ! go directly to the next variable
+   endif
 
 !==============================================================================
@@ -527 +487 @@
 ! 3.3 Write this variable's definition and attributes to the output file
 !==============================================================================
-   if (ndim.ge.3) then
-    units=""
-    title=""
-    ierr=NF_GET_ATT_TEXT(nid,varid,"title",title)
-    ierr=NF_GET_ATT_TEXT(nid,varid,"units",units)
-
-    call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
-   end if
+   units=" "
+   long_name=" "
+   ierr=nf_get_att_text(nid,varid,"long_name",long_name)
+   if (ierr.ne.nf_noerr) then
+   ! if no attribute "long_name", try "title"
+     ierr=nf_get_att_text(nid,varid,"title",long_name)
+   endif
+   ierr=nf_get_att_text(nid,varid,"units",units)
+
+   call def_var(nout,var(j),long_name,units,ndim,dim,varidout,ierr)
 
 !==============================================================================
 ! 3.4 Read variable
-!==== ==========================================================================
-
-   if (ndim.ge.3) then
-      ierr = NF_GET_VAR_REAL(nid,varid,var3d)
-      miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
-      validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
-   end if
-
+!==============================================================================
+
+   ierr = NF_GET_VAR_REAL(nid,varid,var3d)
+   miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
+   validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
 
 !==============================================================================
@@ -578 +537 @@
                endif 
             enddo
-          else   ! average = n
-            call interpolf(timels(n),resultat,lsgcm,var3dxy,timelen)
+          else   ! average = 'n'
+            call interpolf(timels(n),resultat,missing,lsgcm,var3dxy,timelen)
             var3dls(x,y,n,1)=resultat
           endif
@@ -612 +571 @@
                endif 
             enddo
-          else   ! average = n
-             call interpolf(timels(n),resultat,lsgcm,var3dxy,timelen)
+          else   ! average = 'n'
+             call interpolf(timels(n),resultat,missing,lsgcm,var3dxy,timelen)
              var3dls(x,y,k,n)=resultat
-          end if
+          endif
          enddo
         enddo
@@ -626 +585 @@
 !==============================================================================
 
-
-
-   if (ndim.ge.3) then
-
-     ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3dls)
-
-
-     if (ierr.ne.NF_NOERR) then
-      write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
-      stop ""
-     endif
+   ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3dls)
+
+   if (ierr.ne.NF_NOERR) then
+     write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
+     stop ""
+   endif
 
 ! In case there is a "missing value" attribute and "valid range"
-     if (miss.eq.NF_NOERR) then
-      call missing_value(nout,varidout,missing)  
-     end if
-
-     deallocate(var3d)
-     deallocate(var3dls)
-     deallocate(var3dxy)
-   endif !  if (ndim.ge.3) 
+   if (miss.eq.NF_NOERR) then
+     call missing_value(nout,varidout,missing)  
+   endif
+
+   deallocate(var3d)
+   deallocate(var3dls)
+   deallocate(var3dxy)
 
 enddo ! of do j=1,nbvar loop
@@ -666 +619 @@
 !  Writing ctl file to directly read Ls coordinate in grads
 !  (because of bug in grads that refuse to read date like 0089 in .nc files)
-!open(33,file='lslin.ctl')
+write(*,*)""
+write(*,*) "Writing a controle file to directly read Ls coordinate with Grads..."
+
  open(33,file=infile(1:len_trim(infile)-3)//"_Ls.ctl")
  date= (timels(1)-int(timels(1)))*365.
@@ -686 +641 @@
 99 format("TDEF Time ",i5," LINEAR ", i2,a3,i4.4,1x,i2,a2)
     
-   deallocate(timels)
+deallocate(timels)
+   
+write(*,*)""   
+write(*,*) "Program completed !"
+
 contains
 
@@ -876 +835 @@
    ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
 
-   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"title",18,"Control parameters")
+   ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",18,"Control parameters")
 
    ! End netcdf define mode
@@ -1081 +1040 @@
 
 !******************************************************************************
-subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
+subroutine def_var(nid,name,long_name,units,nbdim,dim,nvarid,ierr)
 !==============================================================================
 ! Purpose: Write a variable (i.e: add a variable to a dataset)
-! called "name"; along with its attributes "title", "units"...
+! called "name"; along with its attributes "long_name", "units"...
 ! to a file (following the NetCDF format)
 !==============================================================================
@@ -1102 +1061 @@
 character (len=*), intent(in) :: name
 ! name(): [netcdf] variable's name
-character (len=*), intent(in) :: title
-! title(): [netcdf] variable's "title" attribute
+character (len=*), intent(in) :: long_name
+! long_name(): [netcdf] variable's "long_name" attribute
 character (len=*), intent(in) :: units
 ! unit(): [netcdf] variable's "units" attribute
@@ -1122 +1081 @@
 
 ! Write the attributes
-ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
+ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",len_trim(adjustl(long_name)),adjustl(long_name))
 ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
 
@@ -1182 +1141 @@
 !   print*,'missing_value: valid_range attribution failed'
 !   print*, NF_STRERROR(ierr)
-!   !write(*,*) 'NF_NOERR', NF_NOERR
-!   !CALL abort
 !   stop ""
 !endif
@@ -1193 +1150 @@
    print*, 'missing_value: missing value attribution failed'
    print*, NF_STRERROR(ierr)
-!    WRITE(*,*) 'NF_NOERR', NF_NOERR
-!    CALL abort
    stop ""
 endif
@@ -1296 +1251 @@
 
 end subroutine sol2ls
-!************************************
-
-Subroutine interpolf(x,y,xd,yd,nd)
-!interpolation, give y = f(x) with array xd,yd known, size nd
-! Version with CONSTANT values oustide limits
-! Variable declaration
-! --------------------
-!  Arguments :
-real x,y
-real xd(*),yd(*)
-integer nd
-!  internal
-integer i
-real y_undefined
-                                                                                
-! run
-! ---
-      y_undefined=1.e20
-                                                                                
-      y=0.
-   if ((x.le.xd(1)).and.(x.le.xd(nd))) then
-    if (xd(1).lt.xd(nd)) y = yd(1) !y_undefined
-    if (xd(1).ge.xd(nd)) y = y_undefined ! yd(nd)
-   else if ((x.ge.xd(1)).and.(x.ge.xd(nd))) then
-    if (xd(1).lt.xd(nd)) y =  y_undefined ! yd(nd)
-    if (xd(1).ge.xd(nd)) y = y_undefined ! yd(1)
-    y = yd(nd)
-   else
-    do i=1,nd-1
-     if ( ( (x.ge.xd(i)).and.(x.lt.xd(i+1)) ).or.((x.le.xd(i)).and.(x.gt.xd(i+1))) ) then
-     y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
-     goto 99
-     end if
-    end do
-   end if
-
-!     write (*,*)' x , y' , x,y
-!     do i=1,nd
-!       write (*,*)' i, xd , yd' ,i, xd(i),yd(i)
-!     end do
-!     stop
-                                                                                
-99   continue
-                                                                                
+
+!******************************************************************************
+subroutine interpolf(x,y,missing,xd,yd,nd)
+!==============================================================================
+! Purpose:
+! Yield y=f(x), where xd() end yd() are arrays of known values,
+! using linear interpolation
+! If x is not included in the interval spaned by xd(), then y is set
+! to a default value 'missing'
+! Note:
+! Array xd() should contain ordered (either increasing or decreasing) abscissas
+!==============================================================================
+implicit none
+!==============================================================================
+! Arguments:
+!==============================================================================
+real,intent(in) :: x ! abscissa at which interpolation is to be done
+real,intent(in) :: missing ! missing value (if no interpolation is performed)
+integer :: nd ! size of arrays
+real,dimension(nd),intent(in) :: xd ! array of known absissas
+real,dimension(nd),intent(in) :: yd ! array of correponding values
+
+real,intent(out) :: y ! interpolated value
+!==============================================================================
+! Local variables:
+!==============================================================================
+integer :: i
+
+! default: set y to 'missing'
+y=missing
+
+   do i=1,nd-1
+     if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
+          ((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
+        if ((yd(i).eq.missing).or.(yd(i+1).eq.missing)) then
+          ! cannot perform the interpolation if an encompasing value
+          ! is already set to 'missing'
+        else
+          !linear interpolation based on encompassing values
+          y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
+        endif
+        exit
+     endif
+   enddo
+
+
 end subroutine interpolf
-     
+
+!******************************************************************************
+
 end program lslin