Changeset 1026 for trunk/LMDZ.GENERIC/libf

Timestamp:

Aug 30, 2013, 4:00:31 PM (11 years ago)

Author:

sglmd

Message:

Added a flexible, 2-layer aerosol scenario (initially for Saturn, but can be simply tuned). Called aeroback2lay.

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• aeropacity.F90 (modified) (3 diffs)
• aerosol_mod.F90 (modified) (1 diff)
• callcorrk.F90 (modified) (3 diffs)
• callkeys.h (modified) (2 diffs)
• iniaerosol.F (modified) (2 diffs)
• inifis.F (modified) (1 diff)
• radii_mod.F90 (modified) (2 diffs)
• suaer_corrk.F90 (modified) (1 diff)

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

@@ -67 +67 @@
 
       LOGICAL,SAVE :: firstcall=.true.
-
       REAL CBRT
       EXTERNAL CBRT
@@ -119 +118 @@
         if (iaero_h2so4.ne.0) then
           print*,'iaero_h2so4= ',iaero_h2so4
+        endif
+        if (iaero_back2lay.ne.0) then
+          print*,'iaero_back2lay= ',iaero_back2lay
         endif
 
@@ -356 +358 @@
  
            
+!     ---------------------------------------------------------
+!==================================================================
+!    Two-layer aerosols (unknown composition)
+!    S. Guerlet (2013)
+!==================================================================
+
+      if (iaero_back2lay .ne.0) then
+           iaer=iaero_back2lay
+!       1. Initialization
+            aerosol(1:ngrid,1:nlayermx,iaer)=0.0
+!       2. Opacity calculation
+          DO ig=1,ngrid
+           DO l=1,nlayer-1
+             aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
+             !! 1. below tropospheric layer: no aerosols
+             IF (pplev(ig,l) .gt. pres_bottom_tropo) THEN
+               aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
+             !! 2. tropo layer
+             ELSEIF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
+               aerosol(ig,l,iaer) = obs_tau_col_tropo*aerosol(ig,l,iaer)
+             !! 3. linear transition
+             ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
+               expfactor=log(obs_tau_col_strato/obs_tau_col_tropo)/log(pres_bottom_strato/pres_top_tropo)
+               aerosol(ig,l,iaer)= obs_tau_col_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)*aerosol(ig,l,iaer)
+             !! 4. strato layer
+             ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .gt. pres_top_strato) THEN
+               aerosol(ig,l,iaer)= obs_tau_col_strato*aerosol(ig,l,iaer)
+             !! 5. above strato layer: no aerosols
+             ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN
+               aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
+             ENDIF
+           ENDDO
+          ENDDO
+
+ !       3. Re-normalize to observed total column
+         tau_col(:)=0.0
+         DO l=1,nlayer
+          DO ig=1,ngrid
+               tau_col(ig) = tau_col(ig) &
+                     + aerosol(ig,l,iaer)/(obs_tau_col_tropo+obs_tau_col_strato)
+            ENDDO
+         ENDDO
+
+         DO ig=1,ngrid
+           DO l=1,nlayer-1
+                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
+           ENDDO
+         ENDDO
+
+
+      end if ! if Two-layer aerosols  
+
 
 ! --------------------------------------------------------------------------

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

@@ -13 +13 @@
       integer :: iaero_h2so4 = 0
       logical :: noaero = .false.
+
+! two-layer simple aerosol model
+      integer :: iaero_back2lay = 0
+      REAL :: obs_tau_col_tropo  !! observed total optical depth in the tropospheric layer (visible)
+      REAL :: obs_tau_col_strato !! observed total optical depth in the stratospheric layer (visible)
+      REAL :: pres_bottom_tropo  !! In Pa !   
+      REAL :: pres_top_tropo     !! In Pa
+      REAL :: pres_bottom_strato  !! In Pa
+      REAL :: pres_top_strato     !! In Pa
+      REAL :: size_tropo  !! particle size in the tropospheric layer
+      REAL :: size_strato !! particle size in the stratospheric layer
       
 !==================================================================

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

@@ -14 +14 @@
       use ioipsl_getincom 
       use gases_h
-      use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad
-      use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4
+      use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad,back2lay_reffrad
+      use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4, iaero_back2lay
       USE tracer_h
 
@@ -286 +286 @@
             print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
          endif
-      end do !iaer=1,naerkind
+          if(iaer.eq.iaero_back2lay)then
+            call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
+         endif
+     end do !iaer=1,naerkind
 
 
@@ -302 +305 @@
            reffrad,QREFvis3d,QREFir3d,                             & 
            tau_col,cloudfrac,totcloudfrac,clearsky)                ! get aerosol optical depths
-
+          
 !-----------------------------------------------------------------------
 !     Starting Big Loop over every GCM column

trunk/LMDZ.GENERIC/libf/phystd/callkeys.h

@@ -22 +22 @@
      &   , tracer, mass_redistrib, varactive, varfixed, satval          &
      &   , sedimentation,water,watercond,waterrain                      &
-     &   , aeroco2,aeroh2o,aeroh2so4                                    &
+     &   , aeroco2,aeroh2o,aeroh2so4,aeroback2lay                         &
      &   , aerofixco2,aerofixh2o                                        &
      &   , hydrology, sourceevol, icetstep, albedosnow                  &
@@ -58 +58 @@
       logical sedimentation
       logical water,watercond,waterrain
-      logical aeroco2,aeroh2o,aeroh2so4 
+      logical aeroco2,aeroh2o,aeroh2so4,aeroback2lay
       logical aerofixco2,aerofixh2o
       logical hydrology

trunk/LMDZ.GENERIC/libf/phystd/iniaerosol.F

@@ -12 +12 @@
 c   (CO2 aerosols, dust, water, chemical species, ice...)   
 c
-c   author: Laura Kerber
+c   author: Laura Kerber, S. Guerlet
 c   ------
 c        
@@ -46 +46 @@
       write(*,*) '--- H2SO4 aerosol = ', iaero_h2so4
       
+      if (aeroback2lay) then
+         ia=ia+1
+         iaero_back2lay=ia
+         !! define aerosols parameters (should be in a .def?)
+         !! Saturn: Roman et al. Icarus 2013 2-layer scenario (ISS observations)
+         obs_tau_col_tropo=8.D0       !! observed total optical depth in the tropospheric layer (visible)
+         obs_tau_col_strato=0.08D0    !! observed total optical depth in the stratospheric layer (visible)
+         pres_bottom_tropo= 66000.0   !! In Pa !   
+         pres_top_tropo= 18000.0      !! In Pa
+         pres_bottom_strato= 10000.0  !! In Pa
+         pres_top_strato= 100.0       !! In Pa
+         size_tropo=2.e-6     !! particle size in the tropospheric layer, in meters
+         size_strato=1.e-7    !! particle size in the stratospheric layer
+      endif
+      write(*,*) '--- Two-layer aerosol = ', iaero_back2lay
+
+      
       write(*,*) '=== Number of aerosols= ', ia
       

trunk/LMDZ.GENERIC/libf/phystd/inifis.F

@@ -431 +431 @@
          call getin("aeroh2so4",aeroh2so4)
          write(*,*)" aeroh2so4 = ",aeroh2so4
+         
+         write(*,*)"Radiatively active two-layer aersols?"
+         aeroback2lay=.false.     ! default value
+         call getin("aeroback2lay",aeroback2lay)
+         write(*,*)" aeroback2lay = ",aeroback2lay
+
 
          write(*,*)"Cloud pressure level (with kastprof only):"

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

@@ -76 +76 @@
                nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1 
             endif
-
-            if(iaer.gt.4)then
-               print*,'Error in callcorrk, naerkind is too high (>4).'
+            
+            if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols
+               reffrad(1:ngrid,1:nlayermx,iaer) = 2.e-6
+               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1 
+            endif
+
+
+
+            if(iaer.gt.5)then
+               print*,'Error in callcorrk, naerkind is too high (>5).'
                print*,'The code still needs generalisation to arbitrary'
                print*,'aerosol kinds and number.'
@@ -323 +330 @@
 !==================================================================
 
+!==================================================================
+   subroutine back2lay_reffrad(ngrid,reffrad,nlayer,pplev)
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the effective radii of particles in a 2-layer model
+!
+!     Authors
+!     -------
+!     Sandrine Guerlet (2013)
+!
+!==================================================================
+ 
+      use aerosol_mod   !! Particle sizes and boundaries of aerosol layers defined there
+     Implicit none
+
+#include "callkeys.h"
+#include "dimensions.h"
+#include "dimphys.h"
+
+      integer,intent(in) :: ngrid
+
+      real, intent(out) :: reffrad(ngrid,nlayermx)      ! particle radii 
+      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
+      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
+      REAL :: expfactor
+      INTEGER l,ig
+            
+      reffrad(:,:)=1e-6  !!initialization, not important
+          DO ig=1,ngrid
+            DO l=1,nlayer-1
+              IF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
+                reffrad(ig,l) = size_tropo
+              ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
+                expfactor=log(size_strato/size_tropo) / log(pres_bottom_strato/pres_top_tropo)
+                reffrad(ig,l)= size_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)
+              ELSEIF (pplev(ig,l) .le. pres_bottom_strato) then
+                reffrad(ig,l) = size_strato
+              ENDIF
+            ENDDO
+          ENDDO
+
+   end subroutine back2lay_reffrad
+!==================================================================
+
 
 

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

@@ -167 +167 @@
 ! added by LK
        endif
+
+       if (iaer.eq.iaero_back2lay) then
+         print*, 'naerkind= back2lay', iaer
+
+!     visible
+         file_id(iaer,1) = 'optprop_saturn_vis_n20.dat'
+         lamrefvis(iaer)=0.8E-6  ! 
+!     infrared
+         file_id(iaer,2) = 'optprop_saturn_ir_n20.dat'
+         lamrefir(iaer)=6.E-6  ! 
+! added by SG
+       endif
+       
        
       enddo