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

Timestamp:

Mar 10, 2020, 4:59:49 PM (5 years ago)

Author:

jvatant

Message:

Fix a bug in the 2-layers aerosol scheme in aeropacity.F90
-> It appeared that the normalization for total column optical depth was a

wrong way to go, the actual strato and tropo values weren't matcihing the
input values, and one of the main consequence was to strongly underestimate
the actual strato optical depth compared to input value in most cases.

-> Now we ensure the normalization to the input values for both layers, while

keeping a dependence in dP for each GCM layer.

-> Also add a sanity check ensuring pres_top_tropo > pres_bottom_strato.
--JVO

File:

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

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

@@ -63 +63 @@
 
       real aerosol0, obs_tau_col_aurora, pm, pente_cloud
+      
+      real dp_strato(ngrid)
+      real dp_tropo(ngrid)
 
       INTEGER l,ig,iq,iaer
@@ -371 +374 @@
 !==================================================================
 !    Two-layer aerosols (unknown composition)
-!    S. Guerlet (2013)
+!    S. Guerlet (2013) - Modif by J. Vatant d'Ollone (2020)
 !==================================================================
 
@@ -379 +382 @@
             aerosol(1:ngrid,1:nlayer,iaer)=0.0
 !       2. Opacity calculation
-          DO ig=1,ngrid
-           DO l=1,nlayer-1
+
+
+!       JVO 20 : Modif to have each of the layers (strato and tropo) correctly normalized
+!                Otherwise we previously had the total optical depth correct but for each
+!                separately, so  it didn't match the input values + what's more normalizing
+!                to the sum was making them non-independent : eg changing tau_tropo was
+!                affecting stratopsheric values of optical depth ...
+!
+!                Note that the main consequence of the former version bug was (in most cases)
+!                to strongly underestimate the stratospheric optical depths compared to the
+!                required values, eg, with tau_tropo=10 and tau_strato=0.1, you actually ended
+!                with an actual tau_strato of 1E-4 ... !
+!
+!                NB : Because of the extra transition opacity if the layers are non contiguous,
+!                be aware that at the the bottom we have tau > tau_strato + tau_tropo
+
+         DO ig=1,ngrid
+          dp_tropo(ig)  = 0.D0
+          dp_strato(ig) = 0.D0
+          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)
+               aerosol(ig,l,iaer) = 0.D0
              !! 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)/1.5
-             !! 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)
+               dp_tropo(ig) = dp_tropo(ig) + aerosol(ig,l,iaer)
+             !! 3. linear transition 
+             ! JVO 20 : This interpolation needs to be done AFTER we set strato and tropo (see below)
+             !! 4. strato layer 
+             ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .ge. pres_top_strato) THEN
+               dp_strato(ig) = dp_strato(ig) + 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)
+               aerosol(ig,l,iaer) = 0.D0
              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
+          ENDDO
+         ENDDO
+
+!       3. Re-normalize to the (input) observed (total) column (for each of the layers)
 
          DO ig=1,ngrid
-           DO l=1,nlayer-1
-                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
-           ENDDO
+          DO l=1,nlayer-1
+               IF (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)/dp_tropo(ig)
+               ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
+                 expfactor=log(pplev(ig,l)/pres_top_tropo)/log(pres_bottom_strato/pres_top_tropo)
+                 aerosol(ig,l,iaer) = (obs_tau_col_strato/dp_strato(ig))**expfactor     &
+                                    * (obs_tau_col_tropo/dp_tropo(ig))**(1.0-expfactor) &
+                                    * aerosol(ig,l,iaer)
+               ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .ge. pres_top_strato) THEN
+                 aerosol(ig,l,iaer) = obs_tau_col_strato*aerosol(ig,l,iaer)/dp_strato(ig)
+               ENDIF
+               write(*,*), aerosol(ig,l,iaer)
+            ENDDO
          ENDDO