Changeset 1946 for trunk/LMDZ.TITAN/libf/phytitan

Timestamp:

Jun 13, 2018, 4:18:12 PM (7 years ago)

Author:

jvatant

Message:

+ Fix a major typo in calchim, in the calculation of krate, leading to almost no photodisso ...
+ Also correct kedd profiles to match values from Yelle et al 08 ( slightly overestimated before )
--JVO

File:

• trunk/LMDZ.TITAN/libf/phytitan/calchim.F90 (modified) (5 diffs)

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

@@ -18 +18 @@
   !      
   !          + J. Vatant d'Ollone
-  !               + 02/17 - adaptation for the issu du generic
+  !               + 02/17 - adaptation for the new generic-forked physics
   !               + 01/18 - 03/18 - Major transformations :
   !                   - Upper chemistry fields are now stored in startfi 
   !                     and defined on a pressure grid from Vervack profile
-  !                   - Altitudes above GCM top are computed using correct g(z) 
-  !                     to be coherent with obs and with dz=10km for UV processes
-  !                     but this should be even better if all physiq "under" was doing this !
   !                   - These modifs enables to run chemistry with others resolution than 32x48x55 !
   !                   - Only the actinic fluxes are still read in a 49-lat input but interp. on lat grid
@@ -307 +304 @@
      kedd(:) = 1.e3 ! Default value =/= zero
 
-     ! NB : Eddy coeffs from 1D models (e.g. Lavvas et al 08) in altitude but they're rather linked to pressure
+     ! NB : Eddy coeffs (e.g. Lavvas et al 08, Yelle et al 08) in altitude but they're rather linked to pressure
      !      Below GCM top we have dynamic mixing and for levs < nld=klev-15 the chem. solver ignores diffusion
 
@@ -316 +313 @@
         IF     ( logp.ge.1.0 .and. logp.le.4.0 ) THEN 
               kedd(l) = 10.**(6.0+1.3*(logp-1.0)/3.0)
-     ! 2E7 at 600 km ~ 10-4 mbar
-        ELSEIF ( logp.gt.4.0 .and. logp.le.6.0 ) THEN
-              kedd(l) = 10.**(7.3+0.35*(logp-4.0))
-     ! 1E8 above 900 km ~ 10-6 mbar
-        ELSEIF ( logp.gt.6.0                   ) THEN 
-             kedd(l) = 1.e8
+     ! 2E7 above 600 km ~ 10-4 mbar
+        ELSEIF ( logp.gt.4.0                   ) THEN 
+             kedd(l) = 2.e7
         ENDIF
      ENDDO
@@ -464 +458 @@
            DO i=1,nd_kim+1 ! nd_kim+1 is dissociation of N2 by GCR
 
-              krpddec   =   (   krpd(idec+1,i,ialt  ,klat)   * (1.0-factalt)                   &
-                              + krpd(idec+1,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
-                          * (   krpd(idec+1,i,ialt  ,klat+1) * (1.0-factalt)                   &
-                              + krpd(idec+1,i,ialt+1,klat+1) * factalt       ) * factlat
+                 krpddec   =   (   krpd(idec+1,i,ialt  ,klat)   * (1.0-factalt)                   &
+                                 + krpd(idec+1,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
+                             + (   krpd(idec+1,i,ialt  ,klat+1) * (1.0-factalt)                   &
+                                 + krpd(idec+1,i,ialt+1,klat+1) * factalt       ) * factlat
 
               if      ( factdec.lt.0. ) then 
                  krpddecm1 =   (   krpd(idec  ,i,ialt  ,klat)   * (1.0-factalt)                   &
                                  + krpd(idec  ,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
-                             * (   krpd(idec  ,i,ialt  ,klat+1) * (1.0-factalt)                   &
+                             + (   krpd(idec  ,i,ialt  ,klat+1) * (1.0-factalt)                   &
                                  + krpd(idec  ,i,ialt+1,klat+1) * factalt       ) * factlat
                  krate(l,i) = krpddecm1 * abs(factdec) + krpddec   * ( 1.0 + factdec)
@@ -478 +472 @@
                  krpddecp1 =   (   krpd(idec+2,i,ialt  ,klat)   * (1.0-factalt)                   &
                                  + krpd(idec+2,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
-                             * (   krpd(idec+2,i,ialt  ,klat+1) * (1.0-factalt)                   &
+                             + (   krpd(idec+2,i,ialt  ,klat+1) * (1.0-factalt)                   &
                                  + krpd(idec+2,i,ialt+1,klat+1) * factalt       ) * factlat
                  krate(l,i) = krpddecp1 * factdec      + krpddec   * ( 1.0 - factdec)