Changeset 1663 for trunk/LMDZ.TITAN/libf

Timestamp:

Feb 14, 2017, 11:02:36 AM (8 years ago)

Author:

jvatant

Message:

+ Modifs in haze routine

+ Interpolation corrected (log-lin instead of bilin)
+ Added the case of vertical extension

+ Output useful diagnostics in physiq_mod

Location:

trunk/LMDZ.TITAN/libf/phytitan

Files:

• disr_haze.F90 (modified) (3 diffs)
• physiq_mod.F90 (modified) (1 diff)

trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90

@@ -1 +1 @@
 subroutine disr_haze(dz,press,wno,taeros,ssa,cbar)
+
 use datafile_mod, only: datadir
+
 IMPLICIT NONE
 
-real*8,intent(in)   :: dz,press,wno
-real*8,intent(inout):: taeros,ssa,cbar
+! ==========================================================================
+!
+!  Purpose :
+!    Interpolate values of extinction coefficient, single scattering albedo
+!    and asymetry factor from hazetable ( Lavvas et al. 2010, mean profile, no
+!    detached layer)
+!
+!    + JVO 2017 : Vertical extension out of table implemented
+!
+!  Author :
+!     Jan Vatant d'Ollone (2016)
+!
+! ==========================================================================
+
+real*8,intent(in)   :: dz, press, wno
+real*8,intent(inout):: taeros, ssa, cbar
 
 !---------------------------
-! Attention !!
+! NB !!
 ! taeros is the integrated extinction over the layer
 ! (extinction * thickness of layer)
 !---------------------------
 
-integer             :: i,j,iw,ip,ierr
-real*8              :: wln,factw,factp
-integer,parameter   :: nbwl_PL=328,nblev_PL=162
-real*8,save         :: ext_PL(nblev_PL,nbwl_PL),ssa_PL(nblev_PL,nbwl_PL)
+integer             :: i, j, iw, ip, ierr
+real*8              :: wln, factw, factp
+real*8              :: tmp_p, fact_t, fact_s, fact_c
+integer,parameter   :: nbwl_PL=328, nblev_PL=162
+real*8,save         :: ext_PL(nblev_PL,nbwl_PL), ssa_PL(nblev_PL,nbwl_PL)
 real*8,save         :: asf_PL(nblev_PL,nbwl_PL)
-real*8,save         :: wl_PL(nbwl_PL),press_PL(nblev_PL)
+real*8,save         :: wl_PL(nbwl_PL), press_PL(nblev_PL)
 logical,save        :: firstcall=.true.
 character(len=15)   :: dummy
 
 if (firstcall) then
-   print*,"DISR HAZE - PANAYOTIS LAVVAS"
-
+   print*,"We use DISR haze mean profile from P.Lavvas"
 
 ! read PL table   
@@ -40 +56 @@
    enddo
    close(11)
-   ! convert press_PL into millibar for comparison to press in the generic (wtf press in mbar but ok ... to be modified at some point anyway)
+   ! convert press_PL into millibar for comparison to press in the generic
    press_PL(:)=press_PL(:)*1E-2
+   
+   ! Correction of the table according to Lavvas et al. 2010 (ext coeff = constant < 80 km / 20 mbar ) because of condensation
+!   do i=1,49
+!     ext_PL(i,:) = ext_PL(50,:)
+!   enddo
 
    firstcall=.false.
@@ -65 +86 @@
 enddo
 
+!----------------- Interpolate values from the hazetable  --------------------
+
 if(iw.ne.nbwl_PL) then
-  factw=(wln-wl_PL(iw))/(wl_PL(iw+1)-wl_PL(iw))
+  factw = (wln-wl_PL(iw)) / (wl_PL(iw+1)-wl_PL(iw))
+else
+  factw = 0. ! If we reach the end of the table we keep the last value
 endif
 if(ip.ne.nblev_PL) then
-  factp=(log10(press)-log10(press_PL(ip)))/(log10(press_PL(ip+1))-log10(press_PL(ip)))
+  factp = (press-press_PL(ip)) / (press_PL(ip+1)-press_PL(ip))
 else
-!  factp=(log10(press)-log10(press_PL(nblev_PL)))/(log10(1.e-10)-log10(press_PL(nblev_PL)))
   factp=0.
 endif
 
-if ((iw.ne.nbwl_PL).and.(ip.ne.nblev_PL)) then
-  taeros= ext_PL(ip,iw)  *(1-factw)*(1-factp)+ext_PL(ip+1,iw)  *(1-factw)*factp &
-         +ext_PL(ip,iw+1)*factw    *(1-factp)+ext_PL(ip+1,iw+1)*factw    *factp
-  ssa   = ssa_PL(ip,iw)  *(1-factw)*(1-factp)+ssa_PL(ip+1,iw)  *(1-factw)*factp &
-         +ssa_PL(ip,iw+1)*factw    *(1-factp)+ssa_PL(ip+1,iw+1)*factw    *factp
-  cbar  = asf_PL(ip,iw)  *(1-factw)*(1-factp)+asf_PL(ip+1,iw)  *(1-factw)*factp &
-         +asf_PL(ip,iw+1)*factw    *(1-factp)+asf_PL(ip+1,iw+1)*factw    *factp
-else if ((iw.eq.nbwl_PL).and.(ip.ne.nblev_PL)) then
-  taeros= ext_PL(ip,iw)*(1-factp)+ext_PL(ip+1,iw)*factp 
-  ssa   = ssa_PL(ip,iw)*(1-factp)+ssa_PL(ip+1,iw)*factp 
-  cbar  = asf_PL(ip,iw)*(1-factp)+asf_PL(ip+1,iw)*factp
-else if ((iw.ne.nbwl_PL).and.(ip.eq.nblev_PL)) then
-  taeros= ext_PL(ip,iw)  *(1-factw)*(1-factp) &
-         +ext_PL(ip,iw+1)*factw    *(1-factp)
-  ssa   = ssa_PL(ip,iw)  *(1-factw)*(1-factp) &
-         +ssa_PL(ip,iw+1)*factw    *(1-factp)
-  cbar  = asf_PL(ip,iw)  *(1-factw)*(1-factp) &
-         +asf_PL(ip,iw+1)*factw    *(1-factp)
-else if ((iw.eq.nbwl_PL).and.(ip.eq.nblev_PL)) then
-  taeros= ext_PL(ip,iw)*(1-factp) 
-  ssa   = ssa_PL(ip,iw)*(1-factp) 
-  cbar  = asf_PL(ip,iw)*(1-factp)
+! Lin-Log interpolation : linear on wln, logarithmic on press
+
+taeros = ( ext_PL(ip,iw)*(1.-factw)   + ext_PL(ip,iw+1)  *factw ) ** (1.-factp) &
+        *( ext_PL(ip+1,iw)*(1.-factw) + ext_PL(ip+1,iw+1)*factw ) ** factp
+
+ssa    = ( ssa_PL(ip,iw)*(1.-factw)   + ssa_PL(ip,iw+1)  *factw ) ** (1.-factp) &
+        *( ssa_PL(ip+1,iw)*(1.-factw) + ssa_PL(ip+1,iw+1)*factw ) ** factp
+
+cbar   = ( asf_PL(ip,iw)*(1.-factw)   + asf_PL(ip,iw+1)  *factw ) ** (1.-factp) &
+        *( asf_PL(ip+1,iw)*(1.-factw) + asf_PL(ip+1,iw+1)*factw ) ** factp
+
+! In case of vertical extension over the max of the table
+! We take the scale height on the last 5 levels (more it's not quite log)
+! Arbitray threshold pressure value, just to deal with the last level press=0
+
+if(ip.eq.nblev_PL) then 
+
+  tmp_p = press
+
+  if ( tmp_p .lt. 1.E-15 ) then 
+    tmp_p = 1.E-15
+  endif
+
+  fact_t = log10( ( ext_PL(ip,iw)*(1.-factw) + ext_PL(ip,iw+1)  *factw )   & 
+           / ( ext_PL(ip-5,iw)*(1.-factw) + ext_PL(ip-5,iw+1)  *factw ) )
+  fact_s = log10( ( ssa_PL(ip,iw)*(1.-factw) + ssa_PL(ip,iw+1)  *factw )   & 
+           / ( ssa_PL(ip-5,iw)*(1.-factw) + ssa_PL(ip-5,iw+1)  *factw ) )
+  fact_c = log10( ( asf_PL(ip,iw)*(1.-factw) + asf_PL(ip,iw+1)  *factw )   & 
+           / ( asf_PL(ip-5,iw)*(1.-factw) + asf_PL(ip-5,iw+1)  *factw ) )
+
+  fact_t = fact_t / log10( press_PL(ip) / press_PL(ip-5) )
+  fact_s = fact_s / log10( press_PL(ip) / press_PL(ip-5) )
+  fact_c = fact_c / log10( press_PL(ip) / press_PL(ip-5) )
+
+  taeros = taeros * ( tmp_p / press_PL(ip) ) ** fact_t
+  ssa    = ssa    * ( tmp_p / press_PL(ip) ) ** fact_s
+  cbar   = cbar   * ( tmp_p / press_PL(ip) ) ** fact_c
+
 endif
 

trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

@@ -1278 +1278 @@
       endif ! end of 'enertest'
 
+        ! Temporary inclusions for winds diagnostics.
+        call writediagfi(ngrid,"zdudif","Turbdiff tend. zon. wind","m s-2",3,zdudif)
+        call writediagfi(ngrid,"zdudyn","Dyn. tend. zon. wind","m s-2",3,zdudyn)
 
         ! Temporary inclusions for heating diagnostics.
-        !call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
-        !call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
-        !call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
-        ! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
+        call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
+        call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
+        call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
+        call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
         
         ! For Debugging.
-        !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
+        call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)