Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 2415)
+++ trunk/LMDZ.MARS/README	(revision 2416)
@@ -3164,2 +3164,8 @@
 "tau_pref_gcm" and "tau_pref_scenario" which are repectively the visible
 dust opacity column (at 610 Pa) in the GCM and prescribed by a dust scenario.
+
+== 15/10/2020 == AB
+Correction in solzenangle.F90 : the missing value was not initialized when the attribute was not
+present in the input file.
++ Resolved Ticket #61 : add comments in the code and in the output file to make the program more
+transparent for the users
Index: trunk/LMDZ.MARS/util/solzenangle.F90
===================================================================
--- trunk/LMDZ.MARS/util/solzenangle.F90	(revision 2415)
+++ trunk/LMDZ.MARS/util/solzenangle.F90	(revision 2416)
@@ -26,6 +26,6 @@
 ! filename(): output file name
 ! vartmp(): temporary variable name (read from netcdf input file)
-!character (len=1) :: ccopy
-! ccpy: 'y' or 'n' answer
+character (len=500) :: globatt
+! globatt: [netcdf] global attribute of the output file
 
 character (len=50) :: altlong_name,altunits,altpositive
@@ -100,8 +100,9 @@
 ! var3D_lt(,,,): 4D array to store a field in local time coordinate
 real, dimension(:), allocatable :: lt_gcm
-real, dimension(:), allocatable :: lt_outc
+real, dimension(:), allocatable :: lt_outc ! local lt_out, put in sol decimal value
 real, dimension(:), allocatable :: var_gcm
 real, dimension(:), allocatable :: intsol
 real, dimension(:,:,:), allocatable :: lt_out
+! lt_out: computed local time (hour) corresponding to sza, put in output file
 real :: sza ! solar zenith angle
 character (len=30) :: szastr ! to store the sza value in a string
@@ -237,12 +238,12 @@
 !==============================================================================
 
-      write(*,*) 
-      write(*,*) "opening "//trim(file)//"..."
-      ierr = NF_OPEN(file,NF_NOWRITE,nid)
-      if (ierr.NE.NF_NOERR) then
-         write(*,*) 'ERROR: Pb opening file '//trim(file)
-         write(*,*) NF_STRERROR(ierr)
-         stop ""
-      endif
+write(*,*) 
+write(*,*) "opening "//trim(file)//"..."
+ierr = NF_OPEN(file,NF_NOWRITE,nid)
+if (ierr.NE.NF_NOERR) then
+   write(*,*) 'ERROR: Pb opening file '//trim(file)
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
 
 !==============================================================================
@@ -311,74 +312,74 @@
 
 ! First call; initialize/allocate
-      allocate(lat(latlen),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate lat(latlen)"
-        stop
-      endif
-      allocate(lon(lonlen),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate lon(lonlen)"
-        stop
-      endif
-      allocate(alt(altlen),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate alt(altlen)"
-        stop
-      endif
-      allocate(ctl(ctllen),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate ctl(ctllen)"
-        stop
-      endif
-
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
-#else
-      ierr = NF_GET_VAR_REAL(nid,latvar,lat)
-#endif
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to load latitude"
-        write(*,*) NF_STRERROR(ierr)
-        stop
-      endif
-
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
-#else
-      ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
-#endif
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to load longitude"
-        write(*,*) NF_STRERROR(ierr)
-        stop
-      endif
-
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
-#else
-      ierr = NF_GET_VAR_REAL(nid,altvar,alt)
-#endif
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to load altitude"
-        write(*,*) NF_STRERROR(ierr)
-        stop
-      endif
+   allocate(lat(latlen),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate lat(latlen)"
+     stop
+   endif
+   allocate(lon(lonlen),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate lon(lonlen)"
+     stop
+   endif
+   allocate(alt(altlen),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate alt(altlen)"
+     stop
+   endif
+   allocate(ctl(ctllen),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate ctl(ctllen)"
+     stop
+   endif
+
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
+#else
+   ierr = NF_GET_VAR_REAL(nid,latvar,lat)
+#endif
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to load latitude"
+     write(*,*) NF_STRERROR(ierr)
+     stop
+   endif
+
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
+#else
+   ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
+#endif
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to load longitude"
+     write(*,*) NF_STRERROR(ierr)
+     stop
+   endif
+
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
+#else
+   ierr = NF_GET_VAR_REAL(nid,altvar,alt)
+#endif
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to load altitude"
+     write(*,*) NF_STRERROR(ierr)
+     stop
+   endif
 ! 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)
-
-      if (ctllen .ne. 0) then
-#ifdef NC_DOUBLE
-        ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
-#else
-        ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
-#endif
-        if (ierr.ne.0) then
-          write(*,*) "Error: failed to load controle"
-          write(*,*) NF_STRERROR(ierr)
-          stop
-        endif
-      endif ! of if (ctllen .ne. 0)
+   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)
+
+   if (ctllen .ne. 0) then
+#ifdef NC_DOUBLE
+     ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
+#else
+     ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
+#endif
+     if (ierr.ne.0) then
+       write(*,*) "Error: failed to load controle"
+       write(*,*) NF_STRERROR(ierr)
+       stop
+     endif
+   endif ! of if (ctllen .ne. 0)
 
 !==============================================================================
@@ -430,107 +431,107 @@
    endif
 
-
 !==============================================================================
 ! 2.4.1 Select local times to be saved "Time" in output file
 !==============================================================================
-      write(*,*) "Planet side of the Solar Zenith Angle ?"
-      read(*,*) planetside
-      if ((trim(planetside).ne."morning").and.(trim(planetside).ne."evening")) then
-        write(*,*) "planetside = "//planetside
-        write(*,*) "Error: the planet side must be morning or evening"
-        stop
-      endif
-      write(*,*) "Solar Zenith angle to interpolate ?"
-      write(*,*) "(between 0 and 180 deg, 90deg for the terminator at the surface)"
-      read(*,*) sza
-      if ((sza.lt.0).or.(sza.ge.180)) then
-        write(*,*) "ERROR: the solar zenith angle must lie within [0;180[ deg"
-        stop
-      endif
-      write(*,*) "Beginning date of the simulation file "
-      write(*,*) "(or median date of the simulation for a stats)?"
-      write(*,*) "(i.e. number of sols since Ls=0 at the Time=0.0 in the input file)"
-      read(*,*) starttimeoffset
-
-      if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
-         write(*,*) 'We detect a ""stats"" file'
-         stats=.true.
+   write(*,*) "Planet side of the Solar Zenith Angle ?"
+   read(*,*) planetside
+   if ((trim(planetside).ne."morning").and.(trim(planetside).ne."evening")) then
+     write(*,*) "planetside = "//planetside
+     write(*,*) "Error: the planet side must be morning or evening"
+     stop
+   endif
+   write(*,*) "Solar Zenith angle to interpolate ?"
+   write(*,*) "(between 0 and 180 deg, 90deg for the terminator at the surface)"
+   read(*,*) sza
+   if ((sza.lt.0).or.(sza.ge.180)) then
+     write(*,*) "ERROR: the solar zenith angle must lie within [0;180[ deg"
+     stop
+   endif
+   write(*,*) "Beginning date of the simulation file "
+   write(*,*) "(or median date of the simulation for a stats)?"
+   write(*,*) "(i.e. number of sols since Ls=0 at the Time=0.0 in the input file)"
+   read(*,*) starttimeoffset
+
+   if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
+      write(*,*) 'We detect a ""stats"" file'
+      stats=.true.
 !        We rewrite the time from "stats" from 0 to 1 sol...
-         do it=1,timelen  ! loop temps in
-               time(it) = time(it)/24.
-         end do
-         nsol =1
-      else   ! case of a diagfi or concatnc file 
-        stats=.false.
+      do it=1,timelen  ! loop time input file
+            time(it) = time(it)/24.
+      end do
+      nsol =1
+   else   ! case of a diagfi or concatnc file 
+     stats=.false.
 !       Number of sol in input file
-        nsol = int(time(timelen)) - int(time(1))
-      end if
-      
-      allocate(intsol(nsol),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate intsol(nsol)"
-        stop
-      endif
-      allocate(lt_out(lonlen,latlen,nsol),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate lt_out(lonlen,latlen,nsol)"
-        stop
-      endif
-      do it=1,nsol
-        intsol(it)= int(time(1)) + it-1.
-        do ilon=1,lonlen
-!         For the computation of Ls, it's supposed that the morning/evening
+     nsol = int(time(timelen)) - int(time(1))
+   end if
+
+   allocate(intsol(nsol),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate intsol(nsol)"
+     stop
+   endif
+   allocate(lt_out(lonlen,latlen,nsol),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate lt_out(lonlen,latlen,nsol)"
+     stop
+   endif
+   do it=1,nsol
+     intsol(it)= int(time(1)) + it-1.
+     do ilon=1,lonlen
+!         For the computation of Ls, we try to take into account the rotation time
+!         of Mars it's supposed that the morning/evening
 !         correspond to LT=6h/18h at the given longitude, which is
 !         then reported to a sol value at longitude 0
-          if (trim(planetside).eq."morning") then
-            tmpsol = intsol(it) + (6.-lon(ilon)/15.)/24. + starttimeoffset
-          else !if trim(planetside).eq."evening"
-            tmpsol = intsol(it) + (18.-lon(ilon)/15.)/24. + starttimeoffset            
-          endif
-          call sol2ls(tmpsol,tmpLs)
-          do ilat=1,latlen
+       if (trim(planetside).eq."morning") then
+         tmpsol = intsol(it) + (6.-lon(ilon)/15.)/24. + starttimeoffset
+       else !if trim(planetside).eq."evening"
+         tmpsol = intsol(it) + (18.-lon(ilon)/15.)/24. + starttimeoffset            
+       endif
+       call sol2ls(tmpsol,tmpLs)
+       do ilat=1,latlen
 !           Compute the Local Time of the solar zenith angle at a given longitude
-            call LTsolzenangle(lat(ilat),tmpLs,planetside,sza,miss_lt,tmpLT)
+         call LTsolzenangle(lat(ilat),tmpLs,planetside,sza,miss_lt,tmpLT)
 !           Deduce the sol value at this longitude
-            if (tmpLT.eq.miss_lt) then 
-              ! if there is no LT corresponding to the sza, put a missing value in lt_out
-              lt_out(ilon,ilat,it)=miss_lt
-            else
-              lt_out(ilon,ilat,it)=tmpLT
-            endif
-          enddo
-        enddo
-      enddo
-
-      if (nsol.gt.1) then
-         timelen_tot=timelen
-      else
+         if (tmpLT.eq.miss_lt) then 
+           ! if there is no LT corresponding to the sza, put a missing value in lt_out
+           lt_out(ilon,ilat,it)=miss_lt
+         else
+           lt_out(ilon,ilat,it)=tmpLT
+         endif
+       enddo
+     enddo
+   enddo
+
+   if (nsol.gt.1) then
+      timelen_tot=timelen
+   else
 !        if only 1 sol available, we must add 1 timestep for the interpolation
-         timelen_tot=timelen+1
-      endif      
-      allocate(lt_gcm(timelen_tot),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate lt_gcm(timelen_tot)"
-        stop
-      endif
-      allocate(var_gcm(timelen_tot),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate var_gcm(timelen_tot)"
-        stop
-      endif
-      allocate(lt_outc(nsol),stat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) "Error: failed to allocate lt_outc(nsol)"
-        stop
-      endif
+      timelen_tot=timelen+1
+   endif      
+   allocate(lt_gcm(timelen_tot),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate lt_gcm(timelen_tot)"
+     stop
+   endif
+   allocate(var_gcm(timelen_tot),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate var_gcm(timelen_tot)"
+     stop
+   endif
+   allocate(lt_outc(nsol),stat=ierr)
+   if (ierr.ne.0) then
+     write(*,*) "Error: failed to allocate lt_outc(nsol)"
+     stop
+   endif
 
 !==============================================================================
 ! 2.4.2. Get output file name
 !==============================================================================
-    if (trim(planetside).eq."morning") then 
-      filename=file(1:len_trim(file)-3)//"_MO.nc"
-    else ! if trim(planetside).eq."evening"
-      filename=file(1:len_trim(file)-3)//"_EV.nc"
-    endif
+   if (trim(planetside).eq."morning") then 
+     filename=file(1:len_trim(file)-3)//"_MO.nc"
+   else ! if trim(planetside).eq."evening"
+     filename=file(1:len_trim(file)-3)//"_EV.nc"
+   endif
 
 !==============================================================================
@@ -539,12 +540,12 @@
 
 
-   ! Initialize output file's lat,lon,alt and time dimensions
-    call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
-          altlong_name,altunits,altpositive,&
-          nout,latdimout,londimout,altdimout,timedimout,timevarout)
-
-   ! Initialize output file's aps,bps and phisinit variables
-    call init2(nid,lonlen,latlen,altlen,altdim,&
-               nout,londimout,latdimout,altdimout)
+  ! Initialize output file's lat,lon,alt and time dimensions
+   call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
+         altlong_name,altunits,altpositive,&
+         nout,latdimout,londimout,altdimout,timedimout,timevarout)
+
+  ! Initialize output file's aps,bps and phisinit variables
+   call init2(nid,lonlen,latlen,altlen,altdim,&
+              nout,londimout,latdimout,altdimout)
 
 
@@ -640,9 +641,9 @@
       endif
 
-         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)
+      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)
 
       
@@ -662,64 +663,70 @@
 #endif
 
+      ! If there is no missing value attribute in the input variable, create one
+      if (miss.ne.NF_NOERR) then
+         missing = miss_lt
+         miss = NF_NOERR
+      endif
+
 !==============================================================================
 ! 2.5.4 Interpolation in local time 
 !==============================================================================
 
-        do ilon=1,lonlen
-!          write(*,*) 'processing longitude =', lon(ilon)
-!          Local time at each longitude :
-           do it=1,timelen  ! loop temps in
-               lt_gcm(it) = time(it) + lon(ilon)/360.
-           enddo
-           if (nsol.eq.1) lt_gcm(timelen+1) = lt_gcm(1)+1
-
-           do ilat=1,latlen
+      do ilon=1,lonlen
+!        Recast GCM local time (sol decimal value) at each longitude :
+         do it=1,timelen  ! loop time input file
+             lt_gcm(it) = time(it) + lon(ilon)/360.
+         enddo
+         if (nsol.eq.1) lt_gcm(timelen+1) = lt_gcm(1)+1
+
+         do ilat=1,latlen
+           do it=1,nsol ! loop time local time
+             if (lt_out(ilon,ilat,it).eq.miss_lt) then
+               lt_outc(it)=miss_lt
+             else
+!              Put computed sza local time (hour) into a sol decimal value :
+               lt_outc(it)=intsol(it) + lt_out(ilon,ilat,it)/24.
+             endif
+!            If the computed local time exceeds the GCM input file time bounds,
+!            put a missing value in the local time variable
+             if ((lt_outc(it).gt.lt_gcm(timelen_tot)).OR.(lt_outc(it).lt.lt_gcm(1))) then
+               lt_out(ilon,ilat,it)=miss_lt
+               lt_outc(it)=miss_lt
+             endif
+           enddo ! local time
+
+           if (ndim==3) then
+             do it=1,timelen  ! loop time input file
+                var_gcm(it) = var3d(ilon,ilat,it,1)
+             enddo
+             if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
              do it=1,nsol ! loop time local time
-               if (lt_out(ilon,ilat,it).eq.miss_lt) then
-                 lt_outc(it)=miss_lt
+               if (lt_outc(it).eq.miss_lt) then
+                 var3d_lt(ilon,ilat,it,1)=missing
                else
-                 lt_outc(it)=intsol(it) + lt_out(ilon,ilat,it)/24.
-               endif
-!              If the computed local time exceeds the GCM input file time bounds,
-!              put a missing value in the local time variable
-               if ((lt_outc(it).gt.lt_gcm(timelen_tot)).OR.(lt_outc(it).lt.lt_gcm(1))) then
-                 lt_out(ilon,ilat,it)=miss_lt
-                 lt_outc(it)=miss_lt
+                 call interpolf(lt_outc(it),var3d_lt(ilon,ilat,it,1),&
+                               missing,lt_gcm,var_gcm,timelen_tot)
                endif
              enddo ! local time
 
-             if (ndim==3) then
-               do it=1,timelen  ! loop temps in
-                  var_gcm(it) = var3d(ilon,ilat,it,1)
+           else if  (ndim==4) then
+             do ialt=1,altlen
+               do it=1,timelen ! loop time input file
+                  var_gcm(it) = var3d(ilon,ilat,ialt,it)
                enddo
                if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
                do it=1,nsol ! loop time local time
                  if (lt_outc(it).eq.miss_lt) then
-                   var3d_lt(ilon,ilat,it,1)=missing
+                   var3d_lt(ilon,ilat,ialt,it)=missing
                  else
-                   call interpolf(lt_outc(it),var3d_lt(ilon,ilat,it,1),&
+                   call interpolf(lt_outc(it),var3d_lt(ilon,ilat,ialt,it),&
                                  missing,lt_gcm,var_gcm,timelen_tot)
                  endif
-               enddo ! local time
-
-             else if  (ndim==4) then
-               do ialt=1,altlen
-                 do it=1,timelen ! loop temps in
-                    var_gcm(it) = var3d(ilon,ilat,ialt,it)
-                 enddo
-                 if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
-                 do it=1,nsol ! loop time local time
-                   if (lt_outc(it).eq.miss_lt) then
-                     var3d_lt(ilon,ilat,ialt,it)=missing
-                   else
-                     call interpolf(lt_outc(it),var3d_lt(ilon,ilat,ialt,it),&
-                                   missing,lt_gcm,var_gcm,timelen_tot)
-                   endif
-                 enddo ! local time   
-               enddo ! alt
-             endif
-             
-           enddo ! lat
-        end do ! lon
+               enddo ! local time   
+             enddo ! alt
+           endif
+
+         enddo ! lat
+      end do ! lon
 
 
@@ -739,9 +746,6 @@
       endif
 
-! In case there is a "valid_range" attribute
-      ! Write "valid_range" and "missing_value" attributes in output file
-      if (miss.eq.NF_NOERR) then
-         call missing_value(nout,varidout,validr,miss,valid_range,missing)
-      endif
+      ! Write "missing_value" (and "valid_range" if there is one) attributes in output file
+      call missing_value(nout,varidout,validr,miss,valid_range,missing)
 
       ! free var3d() array
@@ -809,4 +813,16 @@
    ierr=nf_close(nid)
 
+! Write output file's global attribute
+globatt = "GCM simulation that has been interpolated at the prescribed "//trim(planetside)//&
+          "-side solar zenith angle of "//trim(adjustl(szastr))//"deg all around the globe. "//&
+          "Note that the variable LTsolzenangle stands for the Local (at lon,lat) True Solar Time corresponding "//&
+          "to the SZA. You can derive the corresponding sol at lon=0deg by combining LTsolzenangle, longitude and "//&
+          "Time (integer sol value at lon=0deg) variables. For any question on the program, please contact "//&
+          "antoine.bierjon@lmd.jussieu.fr"
+! Switch to netcdf define mode
+ierr = NF_REDEF (nout)
+ierr = NF_PUT_ATT_TEXT(nout,NF_GLOBAL,'comment',len_trim(globatt),trim(globatt))
+! End netcdf define mode
+ierr = NF_ENDDEF(nout)
 ! Close output file
 ierr=NF_CLOSE(nout)
@@ -941,4 +957,10 @@
 call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
              (/timedimout/),timevarout,ierr)
+! Switch to netcdf define mode
+ierr = NF_REDEF (nout)
+ierr = NF_PUT_ATT_TEXT(nout,timevarout,'comment',135,&
+      "The true unit of the variable Time is the integer value of sol at lon=0deg. A false unit is put to enable reading from Grads or Ferret")
+! End netcdf define mode
+ierr = NF_ENDDEF(nout)
 
 !==============================================================================
