Changeset 1161 for trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

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)

File:

• trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 (modified) (8 diffs)

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