Changeset 1409 for trunk/LMDZ.MARS/util/lslin.F90

Timestamp:

Apr 7, 2015, 11:55:58 AM (10 years ago)

Author:

emillour

Message:

Mars GCM:
Update of lslin utility: Added possibility to bin data (instead of interpolating) over the time intervals.
FF

File:

• trunk/LMDZ.MARS/util/lslin.F90 (modified) (22 diffs)

trunk/LMDZ.MARS/util/lslin.F90

@@ -11 +11 @@
 ! from hard-coded values to a computed value) 
 ! Read controle field, if available TN, October 2013
+! Allow to chose a starting Ls and to average within Ls timestep a
+! (great to compare with TES and MCS data. F. Forget december 2013)
 
 
@@ -23 +25 @@
 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=100) :: outfile
 ! outfile(): output file name
@@ -29 +33 @@
 character (len=1) :: answer
 ! answer: the character 'y' (or 'Y'), if input file is of 'concatnc' type
-integer :: nid,varid,ierr
+character (len=1) :: average
+! average: the character 'y' to average within the Ls timestep
+integer :: nid,varid,ierr,miss
 ! nid: [netcdf] ID # of input file
 ! varid: [netcdf] ID # of a variable
@@ -36 +42 @@
 ! nout: [netcdf] ID # of output file
 ! varidout: [netcdf] ID # of a variable (to write in the output file)
-integer :: i,j,k,l,x,y
+integer :: i,j,k,l,x,y,n
 ! counters for various loops
 integer :: start_var
@@ -56 +62 @@
 ! lsgcm(): time coordinate (in unevenly spaced "Ls")
 ! timels(): new time coordinates (evenly spaced "Ls"; written to output file)
-integer :: latlen,lonlen,altlen,timelen,Nls
+integer :: latlen,lonlen,altlen,timelen,Nls,Sls
 ! latlen: # of elements of lat() array
 ! lonlen: # of elements of lon() array
@@ -91 +97 @@
 real, dimension(:), allocatable :: var3dxy
 ! var3dxy(): to store the temporal evolution of a variable (at fixed lat/lon/alt)
-real :: deltatsol,deltals,resultat
+real :: deltatsol,deltals,resultat,ls0
 ! deltatsol: time step (in sols) of input file data
 ! deltals: time step (in Ls) for the data sent to output
+! ls0: first Ls value for the data sent to output
 ! resultat: to temporarily store the result of the interpolation
 character (len=3) :: mon
@@ -99 +106 @@
 real :: date
 ! date: used to compute/build a fake date (for grads output)
+real :: missing
+! to handle missing value when reading /writing files
+real, dimension(2) :: valid_range
+! valid range
 
 !==============================================================================
@@ -229 +240 @@
    else
       ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
+      units_alt="                                                    "
+      ierr=NF_GET_ATT_TEXT(nid,varid,"units",units_alt)
    endif
    write(*,*) "altlen: ",altlen
@@ -256 +269 @@
 
 call initiate(outfile,lat,lon,alt,nout,&
-           latdimout,londimout,altdimout,timedimout,timevarout)
+           latdimout,londimout,altdimout,timedimout,timevarout,units_alt)
 
 !==============================================================================
@@ -334 +347 @@
 do i=1,timelen
    call sol2ls(day_ini+time(i),lsgcm(i)) 
+  if (lsgcm(i).lt.lsgcm(1)) lsgcm(i)=lsgcm(i) + 360.
 enddo
 
@@ -356 +370 @@
     if (0.6*deltatsol.ge.3) deltals=3.
     if (0.6*deltatsol.ge.5) deltals=5.
+    ls0=lsgcm(1) ! First Ls date
 else
     write(*,*) "New timestep in Ls (deg)"
     read(*,*) deltals
-endif
-  NLs=int(Int((lsgcm(timelen)-lsgcm(1))/deltals) +1)
+    write(*,*) "First Ls (deg)"
+    read(*,*) ls0
+    if (ls0.lt.lsgcm(1)) then 
+       write(*,*) 'with this file, the earliest Ls is ',lsgcm(1),'let s use it'
+       ls0=lsgcm(1)
+    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
 !********************************************
 allocate(timels(Nls))
 
 ! Build timels()
-timels(1) = 0.01*nint(100.*lsgcm(1)) ! Ehouarn: what the !!!???!!!
+timels(1) = 0.01*nint(100.*ls0) ! Ehouarn: what the !!!???!!!
 do k=2,Nls
    timels(k) = timels(k-1) + deltals
@@ -470 +497 @@
 
    ierr = NF_GET_VAR_REAL(nid,varid,var3d)
-
-
-!==============================================================================
-! 3.5 Interpolate
+   miss = NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
+   miss = NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
+
+
+!==============================================================================
+! 3.6 interpolate or average
 !==============================================================================
 
@@ -483 +512 @@
 !        write(*,*) 'd: x, y', x, y
         do l=1,timelen
-        var3dxy(l)=var3d(x,y,l,1)
+          var3dxy(l)=var3d(x,y,l,1)
         enddo
-        do l=1,Nls
-        call interpolf(timels(l),resultat,lsgcm,var3dxy,timelen)
-        var3dls(x,y,l,1)=resultat
+        do n=1,Nls
+          if(average.eq.'y') then
+            resultat=0.
+            Sls=0 ! (gcm data counter within each Ls timestep) 
+            do l=1,timelen
+               if((lsgcm(l).ge.(timels(n)-deltals/2.)).and.   &
+               (lsgcm(l).lt.(timels(n)+deltals/2.))) then
+                  if(var3dxy(l) .ne. missing) then
+                    Sls= Sls +1 
+                    resultat = resultat + var3dxy(l)
+                  end if 
+               end if 
+               if (Sls.ne.0) then
+                  var3dls(x,y,n,1)=resultat/float(Sls)
+               else
+                  var3dls(x,y,n,1)=missing
+               endif 
+            enddo
+          else   ! average = n
+            call interpolf(timels(n),resultat,lsgcm,var3dxy,timelen)
+            var3dls(x,y,n,1)=resultat
+          endif
         enddo
        enddo
@@ -500 +548 @@
         var3dxy(l)=var3d(x,y,k,l)
         enddo
-         do l=1,Nls
-        call interpolf(timels(l),resultat,lsgcm,var3dxy,timelen)
-        var3dls(x,y,k,l)=resultat
+        do n=1,Nls
+          if(average.eq.'y') then
+            resultat=0.
+            Sls=0 ! (gcm data counter within each Ls timestep) 
+            do l=1,timelen
+               if((lsgcm(l).ge.(timels(n)-deltals/2.)).and.   &
+               (lsgcm(l).lt.(timels(n)+deltals/2.))) then
+                  if(var3dxy(l) .ne. missing) then
+                    Sls= Sls +1 
+                    resultat = resultat + var3dxy(l)
+                  end if 
+               end if 
+               if (Sls.ne.0) then
+                  var3dls(x,y,k,n)=resultat/float(Sls)
+               else
+                  var3dls(x,y,k,n)=missing
+               endif 
+            enddo
+          else   ! average = n
+             call interpolf(timels(n),resultat,lsgcm,var3dxy,timelen)
+             var3dls(x,y,k,n)=resultat
+          end if
          enddo
         enddo
@@ -510 +577 @@
 
 !==============================================================================
-! 3.4 Write variable to output file
+! 3.7 Write variable to output file
 !==============================================================================
 
@@ -523 +590 @@
       stop ""
    endif
+
+! In case there is a "missing value" attribute and "valid range"
+   if (miss.eq.NF_NOERR) then
+      call missing_value(nout,varidout,valid_range,missing)  
+   end if
 
    deallocate(var3d)
@@ -561 +633 @@
  if(date.ge.305) mon='nov'
  if(date.ge.335) mon='dec'
-write(33,98) outfile
+write(33,98) "^"//outfile
 98 format("DSET ",a)
   write(33,99) Nls, 15,mon, int(timels(1)),nint(deltals*12),'mo'
@@ -571 +643 @@
 !************************************
 subroutine initiate (filename,lat,lon,alt,&
-                     nout,latdimout,londimout,altdimout,timedimout,timevarout)
+                     nout,latdimout,londimout,altdimout,timedimout,timevarout,units_alt)
 !==============================================================================
 ! Purpose:
@@ -607 +679 @@
 integer, intent(out):: timevarout
 ! timevarout: [netcdf] "Time" (considered as a variable) ID
+character (len=50) :: units_alt
+! var(): units of altitude which can be m or Pa (after zrecast)
+
 
 !==============================================================================
@@ -674 +749 @@
 
 ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",8,"altitude")
-ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,"km")
-ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',2,"up")
+!********************* temlporaire ****************
+!ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,units_alt)
+!if((units_alt(1:2).eq.'Pa').or.(units_alt(1:2).eq.'pa')) then
+!   ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',2,"down")
+!else
+!   ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',2,"up")
+!end if
+!********************* temlporaire ****************
+ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,'Pa')
+   ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',4,'down')
+!********************* temlporaire ****************
 
 ! End netcdf define mode
@@ -762 +846 @@
 data year_day /669./            ! # of sols in a martian year
 data peri_day /485.0/           
-data timeperi /1.9082314/
+data timeperi /1.9082314/! True anomaly at vernal equinox = 2*PI-Lsperi
 data e_elips  /0.093358/
 data twopi       /6.2831853/    ! 2.*pi
@@ -872 +956 @@
      
 end
+
+
+
+!******************************************************************************
+!******************************************************************************
+subroutine  missing_value(nout,nvarid,valid_range,missing)
+!==============================================================================
+! Purpose:
+! Write "valid_range" and "missing_value" attributes (of a given
+! variable) to a netcdf file
+!==============================================================================
+! Remarks:
+! NetCDF file must be open
+! Variable (nvarid) ID must be set
+!==============================================================================
+
+implicit none
+
+include "netcdf.inc"  ! NetCDF definitions
+
+!==============================================================================
+! Arguments:
+!==============================================================================
+integer, intent(in) :: nout
+! nout: [netcdf] file ID #
+integer, intent(in) :: nvarid
+! varid: [netcdf] variable ID #
+real, dimension(2), intent(in) :: valid_range
+! valid_range(2): [netcdf] "valid_range" attribute (min and max)
+real, intent(in) :: missing
+! missing: [netcdf] "missing_value" attribute
+
+!==============================================================================
+! Local variables:
+!==============================================================================
+integer :: ierr
+! ierr: [netcdf] subroutine returned error code
+!      INTEGER nout,nvarid,ierr
+!      REAL missing
+!      REAL valid_range(2)
+
+! Switch to netcdf dataset define mode
+ierr = NF_REDEF (nout)
+if (ierr.ne.NF_NOERR) then
+   print*,'missing_value: '
+   print*, NF_STRERROR(ierr)
+endif
+
+! Write valid_range() attribute
+#ifdef NC_DOUBLE
+ierr=NF_PUT_ATT_DOUBLE(nout,nvarid,'valid_range',NF_DOUBLE,2,valid_range)
+#else
+ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
+#endif
+
+if (ierr.ne.NF_NOERR) then
+   print*,'missing_value: valid_range attribution failed'
+   print*, NF_STRERROR(ierr)
+   !write(*,*) 'NF_NOERR', NF_NOERR
+   !CALL abort
+   stop ""
+endif
+
+! Write "missing_value" attribute
+#ifdef NC_DOUBLE
+ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
+#else
+ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
+#endif
+
+if (ierr.NE.NF_NOERR) then
+   print*, 'missing_value: missing value attribution failed'
+   print*, NF_STRERROR(ierr)
+!    WRITE(*,*) 'NF_NOERR', NF_NOERR
+!    CALL abort
+   stop ""
+endif
+
+! End netcdf dataset define mode
+ierr = NF_ENDDEF(nout)
+
+end subroutine  missing_value
+!******************************************************************************