Changeset 3584 for trunk/LMDZ.COMMON/libf/evolution/pem.F90

Timestamp:

Jan 17, 2025, 5:14:18 PM (2 weeks ago)

Author:

jbclement

Message:

PEM:

• Albedo is now updated only at the end of the PEM (and not at every iteration) + Correct way to set it taking into account CO2/H2O ice and frost.
• Cosmetic cleanings.

JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/pem.F90 (modified) (5 diffs)

trunk/LMDZ.COMMON/libf/evolution/pem.F90

@@ -75 +75 @@
 #ifndef CPP_STD
     use comsoil_h,          only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, inertiesoil, flux_geo, nqsoil, qsoil
-    use surfdat_h,          only: tsurf, emis, qsurf, watercap, ini_surfdat_h, &
-                                  albedodat, zmea, zstd, zsig, zgam, zthe,     &
-                                  albedo_h2o_frost,frost_albedo_threshold,     &
-                                  emissiv, watercaptag, perennial_co2ice
+    use surfdat_h,          only: tsurf, qsurf, emis, emissiv, emisice, ini_surfdat_h,   &
+                                  albedodat, albedice, albedo_h2o_frost, albedo_h2o_cap, &
+                                  zmea, zstd, zsig, zgam, zthe, frost_albedo_threshold,  &
+                                  watercap, watercaptag, perennial_co2ice, albedo_perennialco2
     use dimradmars_mod,     only: totcloudfrac, albedo
     use dust_param_mod,     only: tauscaling
@@ -287 +287 @@
 
 ! Loop variables
-integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil
+integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
 
 ! Elapsed time with system clock
@@ -384 +384 @@
     endif
 
-    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,            &
-                         time_0,ps_start_dyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
+    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,                 &
+                         time_0,ps_start_dyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
                          play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
     ps_start_dyn(2) = ps_start_dyn(1)
@@ -844 +844 @@
                         enddo
                     enddo outer2
-                endif
-                if (h2o_ice(ig,islope) > frost_albedo_threshold) then
-                    albedo(ig,1,islope) = albedo_h2o_frost
-                    albedo(ig,2,islope) = albedo_h2o_frost
-                else
-                    albedo(ig,1,islope) = albedodat(ig)
-                    albedo(ig,2,islope) = albedodat(ig)
-                    emis(ig,islope) = emissiv
                 endif
             else if (co2_ice(ig,islope) > 1.e-10 .and. d_co2ice(ig,islope) > 1.e-10) then ! Put tsurf as tcond CO2
@@ -1115 +1107 @@
 deallocate(zplev_start0,zplev_new)
 
+! III_a.6 Albedo update for start file
+do ig = 1,ngrid
+    if (latitude(ig) < 0.) then
+        icap = 2 ! Southern hemisphere
+    else
+        icap = 1 ! Northern hemisphere
+    endif
+    do islope = 1,ngrid
+        ! Bare ground
+        albedo(ig,:,islope) = albedodat(ig)
+        emis(ig,islope) = emissiv
+
+        ! CO2 ice/frost is treated after H20 ice/frost because it is considered dominant
+        ! H2O ice
+        if (h2o_ice(ig,islope) > 0.) then
+            albedo(ig,:,islope) = albedo_h2o_cap
+            emis(ig,islope) = 1.
+        endif
+        ! CO2 ice
+        if (co2_ice(ig,islope) > 0.) then
+            albedo(ig,:,islope) = albedo_perennialco2(icap)
+            emis(ig,islope) = emisice(icap)
+        endif
+        ! H2O frost
+        if (qsurf(ig,igcm_h2o_ice,islope) > 0.) then
+            albedo(ig,:,islope) = albedo_h2o_frost
+            emis(ig,islope) = 1.
+        endif
+        ! CO2 frost
+        if (qsurf(ig,igcm_co2,islope) > 0.) then
+            albedo(ig,:,islope) = albedice(icap)
+            emis(ig,islope) = emisice(icap)
+        endif
+    enddo
+enddo
+
+! III_a.7 Orbital parameters update for start file
 if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)