Changeset 1326

Timestamp:

Aug 19, 2014, 12:42:32 AM (10 years ago)

Author:

jleconte

Message:

18/08/2014 == JL

• corrected insolation calculation in locked and inclined case (proper calculation of sub solar longitude and declination)

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/phystd/comsaison_h.F90 (modified) (1 diff)
• libf/phystd/orbite.F (modified) (6 diffs)
• libf/phystd/physiq.F90 (modified) (5 diffs)

trunk/LMDZ.GENERIC/README

@@ -1045 +1045 @@
 == 18/08/2014 == JL
 - fixed variable type declaration bug in ave_stelspec
+- corrected insolation calculation in locked and inclined case (proper calculation of sub solar longitude and declination)

trunk/LMDZ.GENERIC/libf/phystd/comsaison_h.F90

@@ -6 +6 @@
        integer isaison
        logical callsais
-       real dist_star,declin
+       real dist_star,declin,right_ascen
 
        real, allocatable, dimension(:) :: mu0,fract
-!$OMP THREADPRIVATE(isaison,callsais,dist_star,declin,mu0,fract)
 
        end module comsaison_h

trunk/LMDZ.GENERIC/libf/phystd/orbite.F

@@ -1 @@
-      subroutine orbite(pls,pdist_star,pdecli)
-      
-      use planete_mod, only: p_elips, e_elips, timeperi, obliquit
+      subroutine orbite(pls,pdist_star,pdecli,pright_ascenc)
       implicit none
 
@@ -19 +17 @@
 !     pdist_star    Distance Star-Planet in UA
 !     pdecli        declinaison ( in radians )
+!     pright_ascenc right ascension ( in radians )
 !
 !=======================================================================
@@ -25 +24 @@
 c   -------------
 
-!#include "planete.h"
+#include "planete.h"
 #include "comcstfi.h"
 
@@ -31 +30 @@
 c ----------
 
-      REAL pday,pdist_star,pdecli,pls,i
+      REAL pday,pdist_star,pdecli,pright_ascenc,pls,i
 
 c-----------------------------------------------------------------------
@@ -45 +44 @@
       pdecli = asin (sin(pls)*sin(obliquit*pi/180.))
 
-c ********************* version after 01/01/2000 *******
+c********************* version after 01/01/2000 *******
 c     i=obliquit*pi/180.
 c     pdecli=asin(sin(pls)*sin(i)/sqrt(sin(pls)**2+
@@ -51 +50 @@
 c ******************************************************
 
+c right ascencion
+      If((pls.lt.pi/2.d0)) then
+         pright_ascenc= atan(tan(pls)*cos(obliquit*pi/180.))
+      else if((pls.gt.pi/2.d0).and.(pls.lt.3.d0*pi/2.d0)) then
+         pright_ascenc= pi+atan(tan(pls)*cos(obliquit*pi/180.))
+      else if((pls.gt.3.d0*pi/2.d0)) then
+         pright_ascenc= 2.d0*pi+atan(tan(pls)*cos(obliquit*pi/180.))
+      else if (Abs(pls-pi/2.d0).le.1.d-10) then
+         pright_ascenc= pi/2.d0 
+      else 
+         pright_ascenc=-pi/2.d0 
+      end if
+         
       RETURN
       END

trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

@@ -235 +235 @@
 
       real zls                       ! solar longitude (rad)
+      real zlss                      ! sub solar point longitude (rad)
       real zday                      ! date (time since Ls=0, in martian days)
       real zzlay(ngrid,nlayer)   ! altitude at the middle of the layers
@@ -250 +251 @@
       real zdtdif(ngrid,nlayer)                       ! (K/s)
       real zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer)  ! (m.s-2)
+      real zdvtid(ngrid,nlayer),zdutid(ngrid,nlayer)  ! (m.s-2)
       real zdhadj(ngrid,nlayer)                           ! (K/s)
       real zdtgw(ngrid,nlayer)                            ! (K/s)
@@ -844 +846 @@
 !          Compute local stellar zenith angles
 !          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-           call orbite(zls,dist_star,declin)
+           call orbite(zls,dist_star,declin,right_ascen)
 
            if (tlocked) then
+! JL14 corrects tidally resonant (and inclined) cases. nres=omega_rot/omega_orb
+              zlss=Mod(-(2.*pi*(zday/year_day)*nres - right_ascen),2.*pi)
               ztim1=SIN(declin)
-!              ztim2=COS(declin)*COS(2.*pi*(zday/year_day) - zls*nres)
-!              ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day) - zls*nres)
-! JL12 corrects tidally resonant cases. nres=omega_rot/omega_orb
-              ztim2=COS(declin)*COS(2.*pi*(zday/year_day)*nres - zls)
-              ztim3=-COS(declin)*SIN(2.*pi*(zday/year_day)*nres - zls)
+              ztim2=COS(declin)*COS(zlss)
+              ztim3=COS(declin)*SIN(zlss)
+
 
               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
@@ -858 +860 @@
 
            elseif (diurnal) then
-               ztim1=SIN(declin)
-               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
-               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
-
-               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
+              zlss=-2.*pi*(zday-.5)
+              ztim1=SIN(declin)
+              ztim2=COS(declin)*COS(2.*pi*(zday-.5))
+              ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
+
+              call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
                ztim1,ztim2,ztim3,mu0,fract, flatten)
            else if(diurnal .eqv. .false.) then
@@ -2067 +2070 @@
 
         call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
+        call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
+        call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
+        call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
         call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
         call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)