Changeset 2242 for trunk/LMDZ.TITAN/libf/phytitan/optci.F90

Timestamp:

Feb 17, 2020, 5:44:12 PM (5 years ago)

Author:

jvatant

Message:

+ Update the interface with YAMMS so it now correctly handles the small values of the moments,
requiring dynamics to have a threshold quite low (set to 1D-200) and a sanity check in calmufi
corresponding to this value. Thus it removes 'most' of the unphysical radius obtained in
YAMMS. There are still some, but at least there is no more problem of model stability for the moments
and the code can run.
Still, take care the day you want to calculate opacities from the radii and not the moments.
Although, note that there are some negative values, in the output files, but theses negatives are

harmless, as they are only present in output files, dynamics reseting to epsilon after.

+ Given theses pbs with the radii, also update the optics so that it computes the opacities in
a simpler way, directly for M3, through look-up tables, M3 being a good proxy.
--JVO

File:

• trunk/LMDZ.TITAN/libf/phytitan/optci.F90 (modified) (7 diffs)

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

@@ -6 +6 @@
                          tgasref,pfgasref,wnoi,scalep,indi
   use gases_h
+  use datafile_mod, only: haze_opt_file
   use comcstfi_mod, only: r
   use callkeys_mod, only: continuum,graybody,corrk_recombin,               &
                           callclouds,callmufi,seashaze,uncoupl_optic_haze
   use tracer_h, only : nmicro,nice
-  use MMP_OPTICS
 
   implicit none
@@ -83 +83 @@
   ! variables for k in units m^-1
   real*8 dz(L_LEVELS)
-  !real*8 rho !! see test below
 
   integer igas, jgas, ilay
@@ -89 +88 @@
   integer interm
 
-  real*8 m0as,m3as,m0af,m3af
-  real*8 ext_s,sca_s,ssa_s,asf_s
-  real*8 ext_f,sca_f,ssa_f,asf_f
+  ! Variables for haze optics
+  character(len=200) file_path
+  logical file_ok
+  integer dumch
+  real*8  dumwvl
+
+  real*8 m3as,m3af
+  real*8 dtauaer_s,dtauaer_f
+  real*8,save :: rhoaer_s(L_NSPECTI),ssa_s(L_NSPECTI),asf_s(L_NSPECTI)
+  real*8,save :: rhoaer_f(L_NSPECTI),ssa_f(L_NSPECTI),asf_f(L_NSPECTI)
+!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f)
+
   logical,save :: firstcall=.true.
-  !$OMP THREADPRIVATE(firstcall)
+!$OMP THREADPRIVATE(firstcall)
+
 
   !! AS: to save time in computing continuum (see bilinearbig)
@@ -101 +110 @@
   ENDIF
 
-  ! Some initialisation beacause there's a pb with disr_haze at the limits (nw=1)
+  ! Some initialisation because there's a pb with disr_haze at the limits (nw=1)
   ! I should check this - For now we set vars to zero : better than nans - JVO 2017
-  
-  dhaze_t(:,:) = 0.
-  ssa_t(:,:) = 0.
-  asf_t(:,:) = 0.
+  DHAZE_T(:,:) = 0.0
+  SSA_T(:,:)   = 0.0
+  ASF_T(:,:)   = 0.0
+
+  ! Load tabulated haze optical properties if needed.
+  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
+     OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall
+     READ(12,*) ! dummy header
+     DO NW=1,L_NSPECTI
+       READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw)
+     ENDDO
+     CLOSE(12)
+  ENDIF
 
   !=======================================================================
@@ -118 +137 @@
   do K=2,L_LEVELS
   
-     ilay = k / 2 ! int. arithmetic => gives the gcm layer index
+     ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
      
      DPR(k) = PLEV(K)-PLEV(K-1)
@@ -137 +156 @@
      do K=2,L_LEVELS
      
-        ilay = k / 2 ! int. arithmetic => gives the gcm layer index
+        ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
         
-        ! Optical coupling of YAMMS is plugged but inactivated for now
-        ! as long as the microphysics only isn't fully debugged -- JVO 01/18
         IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
-          m0as = pqmo(ilay,1)
-          m3as = pqmo(ilay,2)
-          m0af = pqmo(ilay,3)
-          m3af = pqmo(ilay,4)
-
-          IF (.NOT.mmp_sph_optics_ir(m0as,m3as,nw,ext_s,sca_s,ssa_s,asf_s)) &
-          CALL abort_gcm("optcv", "Fatal error in mmp_sph_optics_ir", 12)
-          IF (.NOT.mmp_fra_optics_ir(m0af,m3af,nw,ext_f,sca_f,ssa_f,asf_f)) &
-          CALL abort_gcm("optcv", "Fatal error in mmp_fra_optics_ir", 12)
-          dhaze_T(k,nw) = ext_s+ext_f
-          SSA_T(k,nw)   = (sca_s+sca_f)/dhaze_T(k,nw)
-          ASF_T(k,nw)   = (asf_s*sca_s + asf_f*sca_f) /(sca_s+sca_f)
+          m3as = pqmo(ilay,2) / 2.0
+          m3af = pqmo(ilay,4) / 2.0
+          
+          IF ( ilay .lt. 18 ) THEN
+            m3as = pqmo(18,2) / 2.0 *dz(k)/dz(18)
+            m3af = pqmo(18,4) / 2.0 *dz(k)/dz(18)
+          ENDIF
+
+          dtauaer_s     = m3as*rhoaer_s(nw)
+          dtauaer_f     = m3af*rhoaer_f(nw)
+          DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
+          
+          IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN
+            SSA_T(k,nw)   = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
+            ASF_T(k,nw)   = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) )  &
+                            / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f )
+          ELSE
+             DHAZE_T(k,nw) = 0.D0
+             SSA_T(k,nw)   = 1.0
+             ASF_T(k,nw)   = 1.0
+          ENDIF
+          
           IF (callclouds.and.firstcall) & 
             WRITE(*,*) 'WARNING: In optci, optical properties &
@@ -159 +186 @@
         ELSE
           ! Call fixed vertical haze profile of extinction - same for all columns
-          call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
-          if (seashaze) dhaze_T(k,nw) = dhaze_T(k,nw)*seashazefact(k)
+          call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
+          if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
         ENDIF