Changeset 1161 for trunk/LMDZ.GENERIC/libf

Timestamp:

Jan 21, 2014, 4:07:42 PM (11 years ago)

Author:

msylvestre

Message:

LMDZ.GENERIC/UNIVERSAL. rings shadow added in the case diurnal=.false. (daily average of the shadow)

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• callcorrk.F90 (modified) (2 diffs)
• callkeys.h (modified) (2 diffs)
• physiq.F90 (modified) (8 diffs)
• rings.F90 (modified) (6 diffs)

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

@@ -718 +718 @@
 
          if(fract(ig) .ge. 1.0e-4) then ! only during daylight!
-
+!            if(sum(eclipse) .eq. 0.) then
+!                   write(*,*) 'but eclipse is still 0 in physiq !'
+!               endif
             if((ngrid.eq.1).and.(global1d))then
                do nw=1,L_NSPECTV
-                  stel_fract(nw)= stel(nw) * 0.25 / acosz
+                  stel_fract(nw)= stel(nw)* 0.25 / acosz
                                 ! globally averaged = divide by 4
                                 ! but we correct for solar zenith angle
@@ -727 +729 @@
             else
                do nw=1,L_NSPECTV
-                  stel_fract(nw)= stel(nw) * fract(ig)*(1.-eclipse(ig))
+                  stel_fract(nw)= stel(nw) * fract(ig)
                end do
-            endif
-
+            endif 
             call optcv(dtauv,tauv,taucumv,plevrad,                 &
                  qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &

trunk/LMDZ.GENERIC/libf/phystd/callkeys.h

@@ -5 +5 @@
       COMMON/callkeys/callrad,corrk,calldifv,UseTurbDiff,calladj        &
      &   , co2cond,callsoil                                             &
-     &   , season,diurnal,tlocked,iradia,lwrite                         &
+     &   , season,diurnal,tlocked,rings_shadow,iradia,lwrite            &
      &   , iaervar,iddist,topdustref,callstats,calleofdump              &
      &   , enertest                                                     &
@@ -39 +39 @@
      &   , cloudlvl                                                     &
      &   , pceil                                                        &
-     &   , strictboundcorrk                                             &
-     &   , rings_shadow                                         
+     &   , strictboundcorrk                                                                                     
 
       logical callrad,corrk,calldifv,UseTurbDiff                        &
      &   , calladj,co2cond,callsoil                                     &
-     &   , season,diurnal,tlocked,lwrite                                &
+     &   , season,diurnal,tlocked,rings_shadow,lwrite                                &
      &   , callstats,calleofdump                                        &
      &   , callgasvis,continuum,H2Ocont_simple,graybody                 &
-     &   , strictboundcorrk                                             &
-     &   , rings_shadow
+     &   , strictboundcorrk                                             
 
       logical enertest

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

@@ -199 +199 @@
       character*80 fichier 
       integer l,ig,ierr,iq,iaer
-
+      
       !!! this is saved for diagnostic
       real,dimension(:),allocatable,save :: fluxsurf_lw ! incident LW (IR) surface flux (W.m-2)
@@ -411 +411 @@
       logical,save :: ice_update
 
+!     included by MS to compute the daily average of rings shadowing
+      integer, parameter :: nb_hours = 1536 ! set how many times per day are used
+      real :: pas
+      integer :: m
+      ! temporary variables computed at each step of this average 
+      real :: ptime_day ! Universal time in sol fraction 
+      real :: tmp_zls   ! solar longitude
+      real, dimension(:), allocatable :: tmp_fract ! day fraction of the time interval 
+      real, dimension(:), allocatable :: tmp_mu0 ! equivalent solar angle
 !=======================================================================
 
@@ -455 +464 @@
         ALLOCATE(tau_col(ngrid))
 
-         !! this is defined in comsaison_h
-         ALLOCATE(mu0(ngrid))
-         ALLOCATE(fract(ngrid))
-
+        !! this is defined in comsaison_h
+        ALLOCATE(mu0(ngrid))
+        ALLOCATE(fract(ngrid))
+
+            
          !! this is defined in radcommon_h
          ALLOCATE(eclipse(ngrid))
@@ -681 +691 @@
               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)
+!              ztim3=-COS(declin)*SIN(2.  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)
@@ -696 +706 @@
                call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
                ztim1,ztim2,ztim3,mu0,fract)
-
-           else
+           else if(diurnal .eq. .false.) then
 
                call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
-
-           endif
-
+!               write(*,*) 'mu0 =', mu0
+!               write(*,*) 'hello, l 738'
+               ! WARNING: this function appears not to work in 1D
+
+           endif 
+           
            !! Eclipse incoming sunlight (e.g. Saturn ring shadowing)
-           if (rings_shadow) then
-             if (diurnal) then
-               call rings(ngrid, declin, ptime, rad, 0., eclipse)
-               call writediagfi(ngrid,"shad","rings"," ", 2, 1.-eclipse)
-               write(*,*) "Ring shadow activated."
-             else
-               write(*,*) "Ring shadow + diurnal=F not supported yet."
-             endif
-           else
-               eclipse(:) = 0.0
-           endif
+
+           if(rings_shadow) then 
+                write(*,*) 'Rings shadow activated'
+                if(diurnal .eq. .false.) then ! we need to compute the daily average insolation 
+                    pas = 1./nb_hours
+                    ptime_day = 0.
+                    fract(:) = 0.
+                    ALLOCATE(tmp_fract(ngrid))
+                    ALLOCATE(tmp_mu0(ngrid))
+                    tmp_fract(:) = 0.
+                    eclipse(:) = 0.
+                    tmp_mu0(:) = 0.
+                    
+                    do m=1, nb_hours
+                        ptime_day = m*pas
+                        call stellarlong(pday+ptime_day,tmp_zls)
+                        call orbite(tmp_zls,dist_star,declin)
+                        ztim1=SIN(declin)
+                        ztim2=COS(declin)*COS(2.*pi*(pday+ptime_day-.5))
+                        ztim3=-COS(declin)*SIN(2.*pi*(pday+ptime_day-.5))
+                  
+                        call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
+                                 ztim1,ztim2,ztim3,tmp_mu0,tmp_fract)
+                  
+                        call rings(ngrid, declin, ptime_day, rad, 0., eclipse)
+                  
+                        fract(:) = fract(:) + (1.-eclipse(:))*tmp_fract(:) !! fract prend en compte l'ombre des anneaux et l'alternance jour/nuit
+                    enddo
+                
+                    DEALLOCATE(tmp_fract)
+                    DEALLOCATE(tmp_mu0)
+            
+                    fract(:) = fract(:)/nb_hours
+                 
+                 else   
+                    call rings(ngrid, declin, ptime, rad, 0., eclipse)
+                    fract(:) = fract(:) * (1.-eclipse)
+                 endif
+            else
+                 eclipse(:) = 0.0
+            endif
+
 
            if (corrk) then
@@ -1431 +1474 @@
       if(corrk)then
 
+         if(ISNAN(sum(fluxtop_dn))) then 
+            write(*,*) 'err fluxtop_dn'
+         else if(ISNAN(sum(area))) then
+            write(*,*) 'err area'
+         else if(ISNAN(totarea)) then
+            write(*,*) 'area is nan'
+         endif
+
          ISR = SUM(area(:)*fluxtop_dn(:))/totarea
          ASR = SUM(area(:)*fluxabs_sw(:))/totarea
@@ -1717 +1768 @@
            call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
            call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
+           call writediagfi(ngrid,"shad","rings"," ", 2, fract)
 !           call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
 !           call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
@@ -1868 +1920 @@
         !! this is defined in comsaison_h
         IF ( ALLOCATED(mu0)) DEALLOCATE(mu0)
+
         IF ( ALLOCATED(fract)) DEALLOCATE(fract)
+
+
+        !! this is defined in radcommon_h
+        IF ( ALLOCATED(eclipse)) DEALLOCATE(eclipse)
 
         !! this is defined in comsoil_h

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

@@ -40 +40 @@
     REAL, DIMENSION(:), ALLOCATABLE:: x, y, z ! cartesian coordinates of the points on the planet
     REAL :: xs, ys, zs                        ! cartesian coordinates of the points of the subsolar point
-    REAL, DIMENSION(:), ALLOCATABLE :: k      ! parameter (?)
+    REAL, DIMENSION(:), ALLOCATABLE :: k
     REAL, DIMENSION(:), ALLOCATABLE :: N      ! parameter to compute cartesian coordinates on a ellipsoidal planet
     REAL, DIMENSION(:), ALLOCATABLE :: r      ! distance at which the incident ray of sun crosses the equatorial plane
@@ -48 +48 @@
     ! equinox --> no shadow (AS: why is this needed?)
     if(declin .eq. 0.) then
-        eclipse = 0.
+        eclipse(:) = 0.
         return 
     endif 
@@ -63 +63 @@
     ALLOCATE(N(ngrid))
     ALLOCATE(r(ngrid))
-    eclipse = 2000.
+    eclipse(:) = 2000.
 
     ! Size of the rings
@@ -85 +85 @@
 
     ! Convert to cartesian coordinates
-    N = rad/sqrt(1-(e**2)*sinlat**2)
-    x = N*coslat*coslon
-    y = N*coslat*sinlon
-    z = N*(1-e**2)*sinlat
+    N(:) = rad / sqrt(1-(e**2)*sinlat(:)**2)
+    x(:) = N(:)*coslat(:)*coslon(:)
+    y(:) = N(:)*coslat(:)*sinlon(:)
+    z(:) = N(:)*(1-e**2)*sinlat(:)
 
 ! 2) LOCATION OF THE SUBSOLAR POINT 
@@ -95 +95 @@
     ! SG: the minus sign is important! ... otherwise subsolar point adopts a reverse rotation
     phi_S = -(ptime - 0.5)*2.*pi 
-    write(*,*) 'subsol point coords : ', declin*180./pi, phi_S*180./pi
+!    write(*,*) 'subsol point coords : ', declin*180./pi, phi_S*180./pi
 
     ! subsolar latitude is declin (declination of the sun)
@@ -106 +106 @@
 ! 3) WHERE DOES THE INCIDENT RAY OF SUN CROSS THE EQUATORIAL PLAN ?
 
-    k = -z/zs
-    r = (k*xs + x)**2 + (k*ys + y)**2 
-    r = sqrt(r)
+    k(:) = -z(:)/zs
+    r(:) = (k(:)*xs + x(:))**2 + (k(:)*ys + y(:))**2 
+    r(:) = sqrt(r(:))
 
 ! 4) SO WHERE ARE THE SHADOW OF THESE RINGS ?
 
     ! Summer hemisphere is not under the shadow of the rings
-    where(lati*declin .gt. 0.)
-       eclipse = 1000.
+    where(lati(:)*declin .gt. 0.)
+       eclipse(:) = 1000.
     end where
 
     ! No shadow of the rings by night
-    where(x*xs + y*ys + z*zs .lt. 0.)
-       eclipse = 1000.
+    where(x(:)*xs + y(:)*ys + z(:)*zs .lt. 0.)
+       eclipse(:) = 1000.
     end where
 
     ! if the incident rays of sun cross a ring, there is a shadow
     do i=1, nb_A 
-        where(r .ge. A_Rint(i) .and. r .le. A_Rext(i) .and. eclipse .ne. 1000.)
-            eclipse = 1. - exp(-tau_A(i)/cos(declin))
+        where(r(:) .ge. A_Rint(i) .and. r(:) .le. A_Rext(i) .and. eclipse(:) .ne. 1000.)
+            eclipse(:) = 1. - exp(-tau_A(i)/cos(declin))
         end where
     end do 
 
     do i=1, nb_B 
-        where(r .ge. B_Rint(i) .and. r .le. B_Rext(i) .and. eclipse .ne. 1000.)
-            eclipse = 1. - exp(-tau_B(i)/cos(declin))
+        where(r(:) .ge. B_Rint(i) .and. r(:) .le. B_Rext(i) .and. eclipse(:) .ne. 1000.)
+            eclipse(:) = 1. - exp(-tau_B(i)/cos(declin))
         end where
     enddo
     
     do i=1, nb_C 
-        where(r .ge. C_Rint(i) .and. r .le. C_Rext(i) .and. eclipse .ne. 1000.)
-            eclipse = 1. - exp(-tau_C(i)/cos(declin))
+        where(r(:) .ge. C_Rint(i) .and. r(:) .le. C_Rext(i) .and. eclipse(:) .ne. 1000.)
+            eclipse(:) = 1. - exp(-tau_C(i)/cos(declin))
         end where
     enddo
 
     ! At the other places and the excluded ones, eclipse is 0. 
-    where(eclipse .eq. 2000. .or. eclipse .eq. 1000.)
-        eclipse = 0. 
+    where(eclipse(:) .eq. 2000. .or. eclipse(:) .eq. 1000.)
+        eclipse(:) = 0. 
     end where