Changeset 3335 for trunk/LMDZ.GENERIC/libf/phystd

Timestamp:

May 18, 2024, 8:21:24 PM (20 months ago)

Author:

emillour

Message:

Generic PCM:
Add reading/writing of surface albedo in (re)startfi.nc to
improve model restartability. For now only the simpler case
of non-spectral dependent surface albedo is handled.
Turned "surfini" in a module in the process.
Unrelated: added missing delarations in kcm1d so it compiles one again.
EM

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• dyn1d/kcm1d.F90 (modified) (1 diff)
• dyn1d/rcm1d.F (modified) (4 diffs)
• phyetat0_mod.F90 (modified) (5 diffs)
• phyredem.F90 (modified) (4 diffs)
• physiq_mod.F90 (modified) (4 diffs)
• surfini.F (modified) (5 diffs)

trunk/LMDZ.GENERIC/libf/phystd/dyn1d/kcm1d.F90

@@ -106 +106 @@
   LOGICAL :: moderntracdef=.false. ! JVO, YJ : modern traceur.def
 
+  character(len=100) :: dt_file
+  integer :: ios
+  integer :: k
+  
   ! --------------
   ! Initialisation

trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F

@@ -7 +7 @@
       use infotrac, only: nqtot, tname
       use tracer_h, only: noms, is_condensable
+      use radinc_h, only : L_NSPECTV
       use surfdat_h, only: albedodat, phisfi, dryness,
      &                     zmea, zstd, zsig, zgam, zthe,
@@ -94 +95 @@
       integer :: i_h2o_ice=0     ! tracer index of h2o ice
       integer :: i_h2o_vap=0     ! tracer index of h2o vapor
-      REAL emis(1)               ! surface layer
+      REAL emis(1)               ! emissivity of surface
+      real :: albedo(1,L_NSPECTV) ! surface albedo in various spectral bands
       REAL q2(llm+1)             ! Turbulent Kinetic Energy
       REAL zlay(llm)             ! altitude estimee dans les couches (km)
@@ -883 +885 @@
       call getin("albedo",albedodat(1))
       write(*,*) " albedo = ",albedodat(1)
+      ! Initialize surface albedo to that of bare ground
+      albedo(1,:)=albedodat(1)
 
       inertiedat(1,1)=400 ! default value for inertiedat
@@ -954 +958 @@
       call physdem1("startfi.nc",nsoilmx,1,llm,nq,
      &                dtphys,time,
-     &                tsurf,tsoil,emis,q2,qsurf,
+     &                tsurf,tsoil,emis,albedo,q2,qsurf,
      &                cloudfrac,totcloudfrac,hice, 
      &                rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)

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

@@ -7 +7 @@
 subroutine phyetat0 (startphy_file, &
                      ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, &
-                     day_ini,time,tsurf,tsoil, &
-                     emis,q2,qsurf,cloudfrac,totcloudfrac,hice, &
+                     day_ini,time,tsurf,tsoil,emis,albedo, &
+                     q2,qsurf,cloudfrac,totcloudfrac,hice, &
                      rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
 
@@ -18 +18 @@
   use tabfi_mod, only: tabfi
   USE tracer_h, ONLY: noms, igcm_h2o_vap
+  USE radinc_h, ONLY: L_NSPECTV
   USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe
   use iostart, only: nid_start, open_startphy, close_startphy, &
@@ -48 +49 @@
   real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
   real,intent(out) :: emis(ngrid) ! surface emissivity
+  real,intent(out) :: albedo(ngrid,L_NSPECTV) ! albedo of the surface
   real,intent(out) :: q2(ngrid,nlayer+1) ! 
   real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
@@ -129 +131 @@
 
 if (startphy_file) then
-  ! Load bare ground albedo:
+  ! Load bare ground albedo: (will be stored in surfdat_h)
   call get_field("albedodat",albedodat,found)
   if (.not.found) then
@@ -143 +145 @@
 write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
              minval(albedodat), maxval(albedodat)
+
+if (startphy_file) then
+  ! Load surface albedo (for now assume it is spectrally homogeneous)
+  call get_field("albedo",albedo(:,1),found)
+  if (.not.found) then
+    write(*,*) modname//": Failed loading <albedo>"
+    write(*,*) " setting it to bare ground albedo"
+    albedo(1:ngrid,1)=albedodat(1:ngrid)
+  endif
+  ! copy value to all spectral bands
+  do i=2,L_NSPECTV
+    albedo(1:ngrid,i)=albedo(1:ngrid,1)
+  enddo
+else
+  ! If no startfi file, use bare ground value
+  do i=1,L_NSPECTV
+    albedo(1:ngrid,i)=albedodat(1:ngrid)
+  enddo
+endif ! of if (startphy_file)
+write(*,*) "phyetat0: Surface albedo <albedo> range:", &
+             minval(albedo), maxval(albedo)
 
 ! ZMEA

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

@@ -134 +134 @@
 
 subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
-                    phystep,time,tsurf,tsoil,emis,q2,qsurf, &
+                    phystep,time,tsurf,tsoil,emis,albedo,q2,qsurf, &
                     cloudfrac,totcloudfrac,hice, &
                     rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
@@ -141 +141 @@
                       put_var, put_field
   use tracer_h, only: noms
+  USE radinc_h, ONLY: L_NSPECTV
 !  use slab_ice_h, only: noceanmx
   USE ocean_slab_mod, ONLY: nslay
@@ -159 +160 @@
   real,intent(in) :: tsoil(ngrid,nsoil)
   real,intent(in) :: emis(ngrid)
+  real,intent(in) :: albedo(ngrid,L_NSPECTV)
   real,intent(in) :: q2(ngrid,nlay+1)
   real,intent(in) :: qsurf(ngrid,nq)
@@ -189 +191 @@
   call put_field("emis","Surface emissivity",emis)
   
+  ! Surface albedo (assume homegeneous spectral albedo for now)
+  call put_field("albedo","Surface albedo",albedo(:,1))
+  
   ! Planetary Boundary Layer
   call put_field("q2","pbl wind variance",q2)

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

@@ -39 +39 @@
       use time_phylmdz_mod, only: ecritphy, iphysiq, nday
       use phyetat0_mod, only: phyetat0
+      use surfini_mod, only: surfini
       use wstats_mod, only: callstats, wstats, mkstats
       use phyredem, only: physdem0, physdem1
@@ -592 +593 @@
          call phyetat0(startphy_file,                                 &
                        ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
-                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
-                       cloudfrac,totcloudfrac,hice,                   &
+                       day_ini,time_phys,tsurf,tsoil,emis,albedo,     &
+                       q2,qsurf,cloudfrac,totcloudfrac,hice,          &
                        rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
 
@@ -626 +627 @@
 !        Initialize albedo calculation.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         albedo(:,:)=0.0
-         albedo_bareground(:)=0.0
-         albedo_snow_SPECTV(:)=0.0
-         albedo_co2_ice_SPECTV(:)=0.0
-         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
+         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,&
+                      albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
 
 !        Initialize orbital calculation.
@@ -2372 +2370 @@
             call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
                           ptimestep,ztime_fin,                    &
-                          tsurf,tsoil,emis,q2,qsurf_hist,         &
+                          tsurf,tsoil,emis,albedo,q2,qsurf_hist,  &
                           cloudfrac,totcloudfrac,hice,            &
                           rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)

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

@@ +1 @@
+      MODULE surfini_mod
+      
+      IMPLICIT NONE
+      
+      CONTAINS
+
       SUBROUTINE surfini(ngrid,nq,qsurf,albedo,albedo_bareground,
      &                   albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
 
-      USE surfdat_h, only: albedodat
-      USE tracer_h, only: igcm_co2_ice
+      USE surfdat_h, only: albedodat ! bare ground albedo
       use planetwide_mod, only: planetwide_maxval, planetwide_minval
-      use radinc_h, only : L_NSPECTV
+      use radinc_h, only : L_NSPECTV ! number of spectral bands in the visible
       use callkeys_mod, only : albedosnow, albedoco2ice
 
@@ -24 +29 @@
       INTEGER,INTENT(IN) :: ngrid
       INTEGER,INTENT(IN) :: nq
-      REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV)
+      REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV)
       REAL,INTENT(OUT) :: albedo_bareground(ngrid)
       REAL,INTENT(OUT) :: albedo_snow_SPECTV(L_NSPECTV)
@@ -45 +50 @@
       DO ig=1,ngrid
          albedo_bareground(ig)=albedodat(ig)
-         DO nw=1,L_NSPECTV
-            albedo(ig,nw)=albedo_bareground(ig)
-         ENDDO
       ENDDO
       call planetwide_minval(albedo_bareground,min_albedo)
@@ -55 +57 @@
 
 
-      ! Step 3 : We modify the albedo considering some CO2 at the surface. We dont take into account water ice (this is made in hydrol after the first timestep) ...
-      if (igcm_co2_ice.ne.0) then
-         DO ig=1,ngrid
-            IF (qsurf(ig,igcm_co2_ice) .GT. 1.) THEN ! This was changed by MT2015. Condition for ~1mm of CO2 ice deposit.
-               DO nw=1,L_NSPECTV
-                  albedo(ig,nw)=albedo_co2_ice_SPECTV(nw)
-               ENDDO
-            END IF   
-         ENDDO   
-      else
-         write(*,*) "surfini: No CO2 ice tracer on surface  ..."
-         write(*,*) "         and therefore no albedo change."
-      endif     
+      ! Step 3 : Surface albedo already loaded from startfi.nc
+      ! merely report vmin/max values here
+
       call planetwide_minval(albedo,min_albedo)
       call planetwide_maxval(albedo,max_albedo)
@@ -74 +66 @@
 
 
-      END
+      END SUBROUTINE surfini
+      
+      END MODULE surfini_mod