Changeset 2289 for trunk/LMDZ.MARS/util

Timestamp:

Apr 16, 2020, 7:11:13 PM (5 years ago)

Author:

abierjon

Message:

Mars GCM:
Add the program solzenangle.F90 to the LMDZ.MARS/util/ directory, as well as a solzenangle.def
This program actually replaces and seals the dark fate of the former program terminator.F90
It can be used to interpolate GCM files at one given solar zenith angle
(on the morning-side or the evening-side) all around the globe (the terminator corresponding to sza = 90deg)
By Alwont Beebak

Location:

trunk/LMDZ.MARS/util

Files:

• solzenangle.F90 (moved) (moved from trunk/LMDZ.MARS/util/terminator.F90) (14 diffs)
• solzenangle.def (moved) (moved from trunk/LMDZ.MARS/util/terminator.def) (2 diffs)

trunk/LMDZ.MARS/util/solzenangle.F90

@@ -1 @@
-program terminator
+program solzenangle
 
 ! ----------------------------------------------------------------------------
 ! Program to redistribute and interpolate the variables at the local time
-! corresponding to one terminator (morning or evening) everywhere
+! corresponding to one solar zenith angle (morning or evening) everywhere
+! It can especially be used for the terminator (sza=90deg)
 ! input : diagfi.nc  / concat.nc / stats.nc kind of files
 ! authors: A.Bierjon,F.Forget,2020
@@ -33 +34 @@
 ! altpositive(): [netcdf] altitude "positive" attribute
 
-character (len=8) :: sunterm
-! sunterm(): type of terminator ("sunrise" for the morning, "sunset" for the evening)
+character (len=7) :: planetside
+! planetside(): planet side of the solar zenith angle ("morning" or "evening")
 
 integer :: nid,ierr,miss,validr
@@ -103 +104 @@
 real, dimension(:), allocatable :: intsol
 real, dimension(:,:,:), allocatable :: lt_out
+real :: sza ! solar zenith angle
+character (len=30) :: szastr ! to store the sza value in a string
 real :: starttimeoffset ! offset (in sols) of sol 0 in input file, wrt Ls=0
 real :: tmpsol ! to temporarily store a sol value
@@ -113 +116 @@
 real :: miss_lt
 PARAMETER(miss_lt=1E+20)
-! miss_lt: [netcdf] missing value to write for the terminator Local Time in the output file
+! miss_lt: [netcdf] missing value to write for the Local Time in the output file
 real, dimension(2) :: valid_range
 ! valid_range(2): [netcdf] interval in which a value is considered valid
@@ -122 +125 @@
 ! 1.1. Get input file name(s)
 !==============================================================================
+write(*,*) "Welcome in the solar zenith angle interpolator program !"
 write(*,*) 
 write(*,*) "which file do you want to use?  (diagfi... stats...  concat...)"
@@ -430 +434 @@
 ! 2.4.1 Select local times to be saved "Time" in output file
 !==============================================================================
-      write(*,*) 'Type of terminator?'
-      read(*,*) sunterm
-      if ((trim(sunterm).ne."sunrise").and.(trim(sunterm).ne."sunset")) then
-        write(*,*) "sunterm = "//sunterm
-        write(*,*) "Error: terminator type must be sunrise or sunset"
+      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
@@ -470 +481 @@
         do ilon=1,lonlen
 !         For the computation of Ls, it's supposed that the morning/evening
-!         terminator happens at LT=6h/18h at the given longitude, which is
+!         correspond to LT=6h/18h at the given longitude, which is
 !         then reported to a sol value at longitude 0
-          if (trim(sunterm).eq."sunrise") then
+          if (trim(planetside).eq."morning") then
             tmpsol = intsol(it) + (6.-lon(ilon)/15.)/24. + starttimeoffset
-          else !if trim(sunterm).eq."sunset"
+          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 terminator at a given longitude
-            call LTterminator(lat(ilat),tmpLs,sunterm,miss_lt,tmpLT)
+!           Compute the Local Time of the solar zenith angle at a given longitude
+            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 terminator, put a missing value in lt_out
+              ! if there is no LT corresponding to the sza, put a missing value in lt_out
               lt_out(ilon,ilat,it)=miss_lt
             else
@@ -517 +528 @@
 ! 2.4.2. Get output file name
 !==============================================================================
-    if (trim(sunterm).eq."sunrise") then 
-      filename=file(1:len_trim(file)-3)//"_TMO.nc"
-    else ! if trim(sunterm).eq."sunset"
-      filename=file(1:len_trim(file)-3)//"_TEV.nc"
+    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
 
@@ -741 +752 @@
 
 !==============================================================================
-! 2.5.6 Write the terminator Local Time variable in output file
-!==============================================================================
-   write(*,*) "variable LTterminator"
+! 2.5.6 Write the solar zenith angle Local Time variable in output file
+!==============================================================================
+   write(*,*) "variable LTsolzenangle"
    ! build dim(),corner() and edges() arrays
    dim(1)=londimout
@@ -761 +772 @@
    edges(4)=1
 
-   tmpvar="LTterminator"
-   units="hour(LTST)"
-   if (trim(sunterm).eq."sunrise") then
-     title="Morning Terminator Local Time"
-   else !if trim(sunterm).eq."sunset"
-     title="Evening Terminator Local Time"            
+   tmpvar="LTsolzenangle"
+   units="hours"
+   write(szastr,*) sza
+   if (trim(planetside).eq."morning") then
+     title="Morning-side LTST at sza="//trim(adjustl(szastr))//"deg"
+   else !if trim(planetside).eq."evening"
+     title="Evening-side LTST at sza="//trim(adjustl(szastr))//"deg"            
    endif
    call def_var(nout,tmpvar,title,units,3,dim,varidout,ierr)
@@ -802 +814 @@
 write(*,*) "Program completed !"
 
-end program terminator
+end program solzenangle
 
 
@@ -1414 +1426 @@
 
 !******************************************************************************
-subroutine LTterminator(lat,ls,sunterm,missval,lt_term)
-! compute the terminator localtime by inversing
+subroutine LTsolzenangle(lat,ls,planetside,sza,missval,lt_sza)
+! compute the localtime by inversing
 ! the solar zenith angle equation :
-! mu0=sin(lat)*sin(declin)+
-!     cos(lat)*cos(declin)*cos(2.*pi*(localtime/24.-.5))
-! for mu0=cos(sza)=0 (since sza=90deg at the terminator)
+! cos(sza)=sin(lat)*sin(declin)
+!         +cos(lat)*cos(declin)*cos(2.*pi*(localtime/24.-.5))
+
 
 implicit none
@@ -1425 +1437 @@
 real,intent(in) :: lat ! latitude (deg)
 real,intent(in) :: ls ! solar longitude (deg)
-character (len=8), intent(in) :: sunterm ! terminator type : "sunrise" for morning,"sunset" for evening
-real,intent(in) :: missval ! missing value for when there is no terminator (polar day/night)
-real,intent(out) :: lt_term ! local true solar time (hours) of the terminator
+character (len=7), intent(in) :: planetside ! planet side : "morning" or "evening"
+real,intent(in) :: sza ! solar zenith angle (deg) (within [0;180[)
+real,intent(in) :: missval ! missing value for when there is no associated LT (polar day/night)
+real,intent(out) :: lt_sza ! local true solar time (hours) corresponding to the sza
 
 ! local variables:
 double precision :: declin ! declination of the Sun (rad)
+double precision :: tmpcos ! to temporarily store a cosine value
 double precision,parameter :: obliquity=25.1919d0
 double precision,parameter :: pi=3.14159265358979d0
@@ -1439 +1453 @@
 declin=ASIN(sin(ls*degtorad)*sin(obliquity*degtorad))
 
+! Small checks
+if ((cos(lat*degtorad).eq.0).or.(cos(declin).eq.0)) then
+  write(*,*) "ERROR in LTsolzenangle : lat and declin can't have their cosine equal to zero"
+  stop
+endif
+
 ! Compute Local Time
-if (((-tan(lat*degtorad)*tan(declin)).gt.1.).or.(((-tan(lat*degtorad)*tan(declin)).lt.(-1.)))) then
+tmpcos = cos(sza*degtorad)/(cos(lat*degtorad)*cos(declin)) - tan(lat*degtorad)*tan(declin)
+
+if ((tmpcos.gt.1.).or.(tmpcos.lt.(-1.))) then
   ! Detect Polar Day and Polar Night (values out of the arccos definition domain)
-  lt_term = missval
+  lt_sza = missval
   
 else
 
-  if (trim(sunterm).eq."sunrise") then
-  ! Local Time of the Morning Terminator
-    lt_term = 24*(0.5 - ACOS(-tan(lat*degtorad)*tan(declin))/(2.*pi))
-  else ! if trim(sunterm).eq."sunset"
-  ! Local Time of the Evening Terminator
-    lt_term = 24*(0.5 + ACOS(-tan(lat*degtorad)*tan(declin))/(2.*pi))
+  if (trim(planetside).eq."morning") then
+  ! Local Time of the Morning-side Solar Zenith Angle
+    lt_sza = 24*(0.5 - ACOS(tmpcos)/(2.*pi))
+  else ! if trim(planetside).eq."evening"
+  ! Local Time of the Evening-side Solar Zenith Angle
+    lt_sza = 24*(0.5 + ACOS(tmpcos)/(2.*pi))
   endif
   
 endif
 
-end subroutine LTterminator
+end subroutine LTsolzenangle
 
 !******************************************************************************

trunk/LMDZ.MARS/util/solzenangle.def

@@ -4 +4 @@
 temp
 
-sunrise
+morning
+90
 0
 
 -----------------------------------------------------------------------
-ABOVE is the list of inputs to be fed to "terminator.e" if you don t want
+ABOVE is the list of inputs to be fed to "solzenangle.e" if you don t want
 to reply directly to the program:
 
@@ -14 +15 @@
 2) List of X variables to be kept (X lines) or 'all'
 3) blank line at the end
-4) terminator type (sunrise or sunset)
-5) starttimeoffset :
+4) planet side wrt to noon (morning or evening)
+5) solar zenith angle (within [0;180[ deg)
+0deg = zenith ; 90deg = terminator at the surface ; >90deg = terminator in altitude
+6) starttimeoffset :
     - for a diagfi/concat : sol value at the beginning of the run, wrt Ls=0
     - for a stats : sol value at the middle of the run, wrt Ls=0
 
 USE :
-terminator.e < terminator.def
+solzenangle.e < solzenangle.def