Changeset 3131 for trunk/LMDZ.GENERIC

Timestamp:

Nov 20, 2023, 7:52:28 PM (12 months ago)

Author:

yjaziri

Message:

PCM Generic:

Add chemistry escape output + additional flag for photolysis online albedo

YJ

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/aeronostd/calchim_asis.F90 (modified) (4 diffs)
• libf/aeronostd/photolysis_online.F (modified) (2 diffs)

trunk/LMDZ.GENERIC/README

@@ -1828 +1828 @@
 Move mu0 (cosinus of solar zenith angle) calculation out of callrad
 This was a issue for photochemistry because mu0 was updated only iradia times
+
+== 20/11/2023 == YJ
+Add chemistry escape output + additional flag for photolysis online albedo

trunk/LMDZ.GENERIC/libf/aeronostd/calchim_asis.F90

@@ -363 +363 @@
             end do
             if (depos) then
-               ! Surface fluxes
+               ! Surface fluxes and escape
                do iesp=1,nesp
                   call wstats(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux',    & 
@@ -369 +369 @@
                   call wstats(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux',    & 
                                'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
+                  call writediagfi(ngrid,trim(chemnoms(iesp))//'_escape',trim(chemnoms(iesp))//' escape',    &
+                               'molecule.m-2.s-1',2,escape(1,indexchim(trim(chemnoms(iesp)))))
                end do
             endif ! end depos
@@ -395 +397 @@
          end do
          if (depos) then
-            ! Surface fluxes
+            ! Surface fluxes and escape
             do iesp=1,nesp
                call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux_mean',trim(chemnoms(iesp))//' mean flux',    & 
@@ -401 +403 @@
                call writediagfi(ngrid,trim(chemnoms(iesp))//'_flux',trim(chemnoms(iesp))//' flux',    & 
                             'molecule.m-2.s-1',2,surface_flux(1,indexchim(trim(chemnoms(iesp)))))
+               call writediagfi(ngrid,trim(chemnoms(iesp))//'_escape',trim(chemnoms(iesp))//' escape',    &
+                            'molecule.m-2.s-1',2,escape(1,indexchim(trim(chemnoms(iesp)))))
             end do
             ! Restart flux average

trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_online.F

@@ -1172 +1172 @@
 
       use chimiedata_h,  only: albedo_snow_chim, albedo_co2_ice_chim
-!      use slab_ice_h,    only: h_alb_ice, alb_ice_min, alb_ice_max
-      use ocean_slab_mod, only: h_alb_ice
-      use ocean_slab_mod, only: alb_ice_min
-      use ocean_slab_mod, only: alb_ice_max
-      use callkeys_mod,  only: ok_slab_ocean, co2cond, alb_ocean
+      use slab_ice_h,    only: h_alb_ice, 
+     &                         alb_ice_min, alb_ice_max
+      use tracer_h,      only: igcm_h2o_ice, igcm_co2_ice
+      use callkeys_mod,  only: ok_slab_ocean, co2cond, alb_ocean,
+     &                         hydrology
       use phys_state_var_mod, only: albedo_bareground, 
      &                              rnat, qsurf, sea_ice,
@@ -1214 +1214 @@
 !     See hydrol.F90 for taking into account new tendencies
 
-      if(nint(rnat(ig)).eq.0)then
-
-        if(ok_slab_ocean) then
+      if(hydrology)then
+
+        if(nint(rnat(ig)).eq.0)then
+
+          if(ok_slab_ocean) then
         
-          zfra = MAX(0.0,MIN(1.0,qsurf(ig,igcm_h2o_ice)/45.0))     ! Snow Fraction (Critical height 45kg/m2~15cm)
-          alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min)            ! Ice Albedo
-     &                        *exp(-sea_ice(ig)/h_alb_ice)
-          ! Albedo final calculation :
-          do iw=1,nw - 1
-             albedo_chim(iw) = pctsrf_sic(ig)*
-     &                        (albedo_snow_chim(iw)*zfra
-     &                        + alb_ice*(1.0-zfra))
-     &                        + (1.-pctsrf_sic(ig))*alb_ocean
-          enddo  
-
-        else !ok_slab_ocean
+            zfra = MAX(0.0,MIN(1.0,qsurf(ig,igcm_h2o_ice)/45.0))     ! Snow Fraction (Critical height 45kg/m2~15cm)
+            alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min)            ! Ice Albedo
+     &                          *exp(-sea_ice(ig)/h_alb_ice)
+            ! Albedo final calculation :
+            do iw=1,nw - 1
+               albedo_chim(iw) = pctsrf_sic(ig)*
+     &                          (albedo_snow_chim(iw)*zfra
+     &                          + alb_ice*(1.0-zfra))
+     &                          + (1.-pctsrf_sic(ig))*alb_ocean
+            enddo  
+
+          else !ok_slab_ocean
             
             
-!     calculate oceanic ice height including the latent heat of ice formation
-!     hice is the height of oceanic ice with a maximum of maxicethick.
-          hice    = qsurf(ig,igcm_h2o_ice)/rhowater ! update hice to include recent snowfall
-          twater  = tsurf(ig) - hice*RLFTT*rhowater/capcal(ig)
-          ! this is temperature water would have if we melted entire ocean ice layer
-
-          if(twater .lt. T_h2O_ice_liq)then
-
-            do iw=1,nw - 1
-              albedo_chim(iw) = albedo_snow_chim(iw) ! Albedo of ice has been replaced by albedo of snow here. MT2015.
-            enddo
-
+!       calculate oceanic ice height including the latent heat of ice formation
+!       hice is the height of oceanic ice with a maximum of maxicethick.
+            hice    = qsurf(ig,igcm_h2o_ice)/rhowater ! update hice to include recent snowfall
+            twater  = tsurf(ig) - hice*RLFTT*rhowater/capcal(ig)
+            ! this is temperature water would have if we melted entire ocean ice layer
+
+            if(twater .lt. T_h2O_ice_liq)then
+
+              do iw=1,nw - 1
+                albedo_chim(iw) = albedo_snow_chim(iw) ! Albedo of ice has been replaced by albedo of snow here. MT2015.
+              enddo
+
+            else
+
+              DO iw=1,nw - 1
+                albedo_chim(iw) = alb_ocean
+                if(ngrid.eq.1) then
+                  albedo_chim(iw) = albedo_bareground(ig)
+                endif
+              ENDDO
+
+            endif
+
+          endif!(ok_slab_ocean)
+
+
+!       Continent
+!       ---------
+        elseif (nint(rnat(ig)).eq.1) then
+
+!       re-calculate continental albedo
+          if (qsurf(ig,igcm_h2o_ice).ge.snowlayer) then
+            DO iw=1,nw - 1
+              albedo_chim(iw) = albedo_snow_chim(iw)
+            ENDDO
           else
-
             DO iw=1,nw - 1
-              albedo_chim(iw) = alb_ocean
-              if(ngrid.eq.1) then
-                albedo_chim(iw) = albedo_bareground(ig)
-              endif
+              albedo_chim(iw) = albedo_bareground(ig)
+     &                        + (albedo_snow_chim(iw)
+     &                           - albedo_bareground(ig))
+     &                        *qsurf(ig,igcm_h2o_ice)/snowlayer
             ENDDO
-
           endif
 
-        endif!(ok_slab_ocean)
-
-
-!     Continent
-!     ---------
-      elseif (nint(rnat(ig)).eq.1) then
-
-!     re-calculate continental albedo
-        if (qsurf(ig,igcm_h2o_ice).ge.snowlayer) then
-          DO iw=1,nw - 1
-            albedo_chim(iw) = albedo_snow_chim(iw)
-          ENDDO
         else
-          DO iw=1,nw - 1
-            albedo_chim(iw) = albedo_bareground(ig)
-     &                      + (albedo_snow_chim(iw)
-     &                         - albedo_bareground(ig))
-     &                      *qsurf(ig,igcm_h2o_ice)/snowlayer
-          ENDDO
+
+          print*,'Surface type not recognised in photolysis_online.F!'
+          print*,'Exiting...'
+          call abort
+
         endif
-
       else
-
-        print*,'Surface type not recognised in photolysis_online.F!'
-        print*,'Exiting...'
-        call abort
-
-      endif
+        albedo_chim(:) = albedo_bareground(ig)
+      endif ! end if hydrology
 
 !     Re-add the albedo effects of CO2 ice if necessary