Changeset 2297 for trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90

Timestamp:

Apr 26, 2020, 7:45:46 PM (5 years ago)

Author:

jvatant

Message:

Add a generic n-layer aerosol scheme to replace the former buggy 2-layer scheme as well as the hard-coded NH3 cloud.

It can be called using 'aeronlay=.true.' in callphys.def, and set the number of layers (up to 4) with 'nlayaero'.
Then, the following parameters are read as arrays of size nlayaero in callphys.def (separated by blank space)

*aeronlay_tauref (Optical depth of aerosol layer at ref wavelenght)
*aeronlay_lamref (Ref wavelenght (m))
*aeronlay_choice (Choice for vertical profile - 1:tau follows atm scale height btwn top and bottom - 2:tau follows it own scale height)
*aeronlay_pbot (Bottom pressure (Pa))
*aeronlay_ptop (Top pressure (Pa) - useful only if choice=1)
*aeronlay_sclhght (Ratio of aerosol layer scale height / atmospheric scale height - useful only if choice=2 )
*aeronlay_size (Particle size (m))
*optprop_aeronlay_vis File for VIS opt properties.
*optprop_aeronlay_ir File for IR opt properties.

+Extra info :

+ In addition of solving the bug from 2-layer it enables different optical properties.
+ The former scheme are left for retrocompatibility (for now) but you should use the new one.
+ See aeropacity.F90 for the calculations

+ Each layer can have different optical properties, size of particle ...
+ You have different choices for vertical profile of the aerosol layers :

• aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.
• aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).

In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.

+ The reference wavelenght for input optical depth is now read as input (aeronlay_lamref)
+ Following the last point some comment is added in suaer_corrk about the 'not-really-dummy'ness of IR lamref..

File:

• trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90 (modified) (13 diffs)

trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90

@@ -11 +11 @@
                 pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo,     &
                 pres_bottom_strato,pres_top_strato,obs_tau_col_strato,  &
-                tau_nh3_cloud, pres_nh3_cloud 
+                tau_nh3_cloud, pres_nh3_cloud,                          &
+                nlayaero, aeronlay_tauref, aeronlay_choice,             & 
+                aeronlay_pbot, aeronlay_ptop, aeronlay_sclhght
                   
        implicit none
@@ -27 +29 @@
 !     update J.-B. Madeleine (2008)
 !     dust removal, simplification by Robin Wordsworth (2009)
+!     Generic n-layer aerosol - J. Vatant d'Ollone (2020)
 !
 !     Input
@@ -62 +65 @@
       logical,intent(in) :: clearsky
 
-      real aerosol0, obs_tau_col_aurora, pm, pente_cloud
-      
+      real aerosol0, obs_tau_col_aurora, pm
+      real pcloud_deck, cloud_slope
+
       real dp_strato(ngrid)
       real dp_tropo(ngrid)
-
-      INTEGER l,ig,iq,iaer
+      real dp_layer(ngrid)
+
+      INTEGER l,ig,iq,iaer,ia
 
       LOGICAL,SAVE :: firstcall=.true.
@@ -129 +134 @@
         if (iaero_nh3.ne.0) then
           print*,'iaero_nh3= ',iaero_nh3
+        endif
+        if (iaero_nlay(1).ne.0) then
+          print*,'iaero_nlay= ',iaero_nlay(:)
         endif
         if (iaero_aurora.ne.0) then
@@ -375 +383 @@
 !    Two-layer aerosols (unknown composition)
 !    S. Guerlet (2013) - Modif by J. Vatant d'Ollone (2020)
+!    
+!    This scheme is deprecated and left for retrocompatibility
+!    You should use the n-layer scheme below !
+!
 !==================================================================
 
@@ -435 +447 @@
                  aerosol(ig,l,iaer) = obs_tau_col_strato*aerosol(ig,l,iaer)/dp_strato(ig)
                ENDIF
-               write(*,*), aerosol(ig,l,iaer)
             ENDDO
          ENDDO
@@ -445 +456 @@
 !    Saturn/Jupiter ammonia cloud = thin cloud (scale height 0.2 hard coded...)
 !    S. Guerlet (2013)
+!    JVO 20 : You should now use the generic n-layer scheme below
 !==================================================================
 
@@ -461 +473 @@
 
              ELSEIF (pplev(ig,l) .le. pres_nh3_cloud ) THEN 
-             pente_cloud=5. !!(hard-coded, correspond to scale height 0.2)
-             aerosol(ig,l,iaer) = ((pplev(ig,l)/pres_nh3_cloud)**(pente_cloud))*tau_nh3_cloud 
+             cloud_slope=5. !!(hard-coded, correspond to scale height 0.2)
+             aerosol(ig,l,iaer) = ((pplev(ig,l)/pres_nh3_cloud)**(cloud_slope))*tau_nh3_cloud 
 
              ENDIF
@@ -470 +482 @@
           
 !       3. Re-normalize to observed total column
-         tau_col(:)=0.0
+         dp_layer(:)=0.0
          DO l=1,nlayer
           DO ig=1,ngrid
-               tau_col(ig) = tau_col(ig) &
+               dp_layer(ig) = dp_layer(ig) &
                      + aerosol(ig,l,iaer)/tau_nh3_cloud
             ENDDO
@@ -480 +492 @@
          DO ig=1,ngrid
            DO l=1,nlayer-1
-                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
+                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/dp_layer(ig)
            ENDDO
          ENDDO
 
      end if ! if NH3 cloud  
+
+!=========================================================================================================
+!    Generic N-layers aerosols/clouds
+!    Author : J. Vatant d'Ollone (2020)
+!    
+!    Purpose: Replaces and extents the former buggy 2-layer scheme as well as hard-coded NH3 cloud
+!    
+!    + Each layer can have different optical properties, size of particle ...
+!    + Enables up to n=4 layers as we apparently cannot run with more scatterers (could be worth checking...)
+!    + You have different choices for vertical profile of the aerosol layers :
+!           * aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.
+!           * aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).
+!                                   In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.
+!           * aeronlay_choice = ... feel free to add more cases  !
+!    + Layers can overlap if needed (if you want a 'transition layer' as in the 2-scheme, just add it)
+!
+!=========================================================================================================
+
+      if (iaero_nlay(1) .ne.0) then
+
+        DO ia=1,nlayaero
+           iaer=iaero_nlay(ia)
+
+!          a. Initialization
+           aerosol(1:ngrid,1:nlayer,iaer)=0.D0
+
+!          b. Opacity calculation
+           
+           ! Case 1 : Follows atmospheric scale height between boundaries pressures
+           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           IF (aeronlay_choice(ia).eq.1) THEN
+
+             dp_layer(:)=0.D0
+             DO ig=1,ngrid
+               DO l=1,nlayer-1
+                 !! i. Opacity follows scale height
+                 IF ( pplev(ig,l).le.aeronlay_pbot(ia)   .AND.          &
+                      pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN
+                   aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
+                   dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)
+                 !! ii. Outside aerosol layer boundaries: no aerosols
+                 ELSE
+                   aerosol(ig,l,iaer) = 0.D0
+                 ENDIF
+               ENDDO
+             ENDDO
+             ! iii. Re-normalize to required total opacity
+             DO ig=1,ngrid
+               DO l=1,nlayer-1
+                 IF ( pplev(ig,l).le.aeronlay_pbot(ia)   .AND.          &
+                      pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN
+                  aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &
+                                     * aeronlay_tauref(ia)
+                 ENDIF
+               ENDDO
+             ENDDO
+
+           ! Case 2 : Follows input scale height
+           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           ELSE IF (aeronlay_choice(ia).eq.2) THEN
+           
+             cloud_slope  = 1.D0/aeronlay_sclhght(ia)
+             pcloud_deck  = aeronlay_pbot(ia)
+             dp_layer(:)  = 0.D0
+
+             DO ig=1,ngrid
+               DO l=1,nlayer-1
+                 !! i. Below cloud layer: no opacity
+                 IF (pplev(ig,l) .gt. pcloud_deck) THEN
+                   aerosol(ig,l,iaer) = 0.D0            
+                 !! ii. Follows scale height above cloud deck
+                 ELSEIF (pplev(ig,l) .le. pcloud_deck) THEN 
+                   aerosol(ig,l,iaer) = ((pplev(ig,l)/pcloud_deck)**(cloud_slope))
+                   dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)
+                 ENDIF
+               ENDDO
+             ENDDO
+             ! iii. Re-normalize to required total opacity
+             DO ig=1,ngrid
+               DO l=1,nlayer-1
+                 IF (pplev(ig,l) .le. pcloud_deck) THEN 
+                  aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &
+                                     * aeronlay_tauref(ia)
+                 ENDIF
+               ENDDO
+             ENDDO
+
+           ENDIF ! aeronlay_choice
+
+          ENDDO ! loop on n aerosol layers
+
+      end if ! if N-layer aerosols
+  
 !==================================================================
 !    Jovian auroral aerosols (unknown composition) NON-GENERIC: vertical and meridional profile tuned to observations
@@ -506 +611 @@
          
  !       3. Meridional distribution, and re-normalize to observed total column
-         tau_col(:)=0.D0
+         dp_layer(:)=0.D0
          DO ig=1,ngrid
           !!Jupiter
@@ -528 +633 @@
 
           DO l=1,nlayer  
-               tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora
+               dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora
           ENDDO
          ENDDO
@@ -534 +639 @@
          DO ig=1,ngrid
            DO l=1,nlayer-1
-                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
+                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/dp_layer(ig)
            ENDDO
          ENDDO