Changeset 1681 for trunk/LMDZ.TITAN/libf

Timestamp:

Mar 30, 2017, 4:54:02 PM (8 years ago)

Author:

jvatant

Message:

Update deftank + correct a bug in extrapolation in haze opacity table routine

File:

• trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90 (modified) (4 diffs)

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

@@ -30 +30 @@
 integer             :: i, j, iw, ip, ierr
 real*8              :: wln, factw, factp
-real*8              :: tmp_p, fact_t, fact_s, fact_c
+real*8              :: tmp_p, fact_t
 integer,parameter   :: nbwl_PL=328, nblev_PL=162
 real*8,save         :: ext_PL(nblev_PL,nbwl_PL), ssa_PL(nblev_PL,nbwl_PL)
@@ -82 +82 @@
 
 !----------------- Interpolate values from the hazetable  --------------------
+if (iw.ne. nbwl_PL) then
+  factw = (wln-wl_PL(iw)) / (wl_PL(iw+1)-wl_PL(iw))
+endif
 
-if(iw.ne.nbwl_PL) then
-  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
+if (ip .ne. nblev_PL) then
   factp = (press-press_PL(ip)) / (press_PL(ip+1)-press_PL(ip))
-else
-  factp=0.
 endif
 
 ! Lin-Log interpolation : linear on wln, logarithmic on press
 
-if ((iw.ne.nbwl_PL).and.(ip.ne.nblev_PL)) then
+if((ip.ne.nblev_PL) .and. (iw.ne.nbwl_PL)) then 
+
 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
@@ -105 +102 @@
 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
-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
+else if ((ip.ne.nblev_PL) .and. (iw.eq.nbwl_PL)) then
 
-cbar   =  asf_PL(ip,iw) ** (1.-factp) * asf_PL(ip+1,iw) ** factp
+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
+
 
 ! 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
+! We do not touch to ssa and cbar and let them at the value of last level
+! (extrap would lead to too dark aerosols)
 
 else if(ip.eq.nblev_PL) then 
@@ -124 +126 @@
   endif
 
-  if (iw.ne.nbwl_PL) then
-    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 ) )
-  else
-    fact_t = log10( ext_PL(ip,iw) / ext_PL(ip-5,iw) )
-    fact_s = log10( ssa_PL(ip,iw) / ssa_PL(ip-5,iw) )
-    fact_c = log10( asf_PL(ip,iw) / asf_PL(ip-5,iw) )
-  endif            
+  if(iw.ne.nbwl_PL) then  
+
+     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 ) )
+
+  else if (iw.eq.nbwl_PL) then
+
+     fact_t = log10( ext_PL(ip,iw) / ext_PL(ip-5,iw) )  
+
+  endif
 
   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