Index: /trunk/LMDZ.GENERIC/changelog.txt =================================================================== --- /trunk/LMDZ.GENERIC/changelog.txt (revision 3334) +++ /trunk/LMDZ.GENERIC/changelog.txt (revision 3335) @@ -1930,2 +1930,9 @@ not 1 this extra computation of radiative tendencies will break model restartability (the famous "1+1=2" requirement). + +== 18/05/2024 == EM +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. Index: /trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/newstart.F =================================================================== --- /trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/newstart.F (revision 3334) +++ /trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/newstart.F (revision 3335) @@ -21,4 +21,5 @@ USE tracer_h, ONLY: igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat + USE radinc_h, only : L_NSPECTV ! number of spectral bands in the visible USE surfdat_h, ONLY: phisfi, albedodat, & zmea, zstd, zsig, zgam, zthe @@ -112,5 +113,7 @@ REAL q2(ngridmx,llm+1) ! REAL rnaturfi(ngridmx) - real alb(iip1,jjp1),albfi(ngridmx) ! albedos + real alb(iip1,jjp1),albfi(ngridmx) ! bare ground albedos + real albedodyn(iip1,jjp1),albedofi(ngridmx) ! surface albedos + real spectral_albedofi(ngridmx,L_NSPECTV) ! spectral surface albedo real,ALLOCATABLE :: ith(:,:,:),ithfi(:,:) ! thermal inertia (3D) real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D) @@ -361,12 +364,13 @@ fichnom = 'startfi.nc' CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx, - . nqtot,day_ini,time, - . tsurf,tsoil,emis,q2,qsurf, !) ! temporary modif by RDW - . cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice, - . sea_ice) + & nqtot,day_ini,time, + & tsurf,tsoil,emis,spectral_albedofi,q2,qsurf, + & cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice, + & sea_ice) ! copy albedo and soil thermal inertia on (local) physics grid do i=1,ngridmx albfi(i) = albedodat(i) + albedofi(i)= spectral_albedofi(i,1) ! assume same albedo at all wavelenghts do j=1,nsoilmx ithfi(i,j) = inertiedat(i,j) @@ -380,4 +384,5 @@ ! to correctly recast things on physics grid) call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,albedofi,albedodyn) call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith) @@ -1662,5 +1667,5 @@ call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot, & dtphys,real(day_ini), - & tsurf,tsoil,emis,q2,qsurf, + & tsurf,tsoil,emis,spectral_albedofi,q2,qsurf, & cloudfrac,totalfrac,hice, & rnat,pctsrf_sic,tslab,tsea_ice,sea_ice) Index: /trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/start2archive.F =================================================================== --- /trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/start2archive.F (revision 3334) +++ /trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/start2archive.F (revision 3335) @@ -21,5 +21,5 @@ use infotrac, only: infotrac_init, nqtot, tname USE comsoil_h - + USE radinc_h, only : L_NSPECTV ! number of spectral bands in the visible ! use slab_ice_h, only: noceanmx USE ocean_slab_mod, ONLY: nslay @@ -73,4 +73,5 @@ REAL,ALLOCATABLE :: qsurf(:,:) REAL emis(ngridmx) + REAL :: albedo(ngridmx,L_NSPECTV) ! spectral surface albedo INTEGER start,length PARAMETER (length = 100) @@ -99,4 +100,5 @@ REAL,ALLOCATABLE :: qsurfS(:,:) REAL emisS(ip1jmp1) + REAL :: albedoS(ngridmx) ! surface albedo assumed same at all wavelengths ! added by FF for cloud fraction setup @@ -245,5 +247,5 @@ CALL phyetat0(.true.,ngridmx,llm,fichnom,0,Lmodif,nsoilmx,nqtot, . day_ini_fi,timefi, - . tsurf,tsoil,emis,q2,qsurf, + . tsurf,tsoil,emis,albedo,q2,qsurf, ! change FF 05/2011 . cloudfrac,totalcloudfrac,hice, @@ -351,4 +353,5 @@ call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,inertiedat,ithS) call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis,emisS) + call gr_fi_dyn(1,ngridmx,iip1,jjp1,albedo(1,1),albedoS) call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S) call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS) @@ -455,4 +458,6 @@ ! & 'kg/m2',2,co2iceS) call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS) + call write_archive(nid,ntime,'albedo','surface albedo',' ', + & 2,albedoS) call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps) call write_archive(nid,ntime,'tsurf','surf T','K',2,tsurfS) Index: /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/kcm1d.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/kcm1d.F90 (revision 3334) +++ /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/kcm1d.F90 (revision 3335) @@ -106,4 +106,8 @@ LOGICAL :: moderntracdef=.false. ! JVO, YJ : modern traceur.def + character(len=100) :: dt_file + integer :: ios + integer :: k + ! -------------- ! Initialisation Index: /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F (revision 3334) +++ /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F (revision 3335) @@ -7,4 +7,5 @@ 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,5 +95,6 @@ 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,4 +885,6 @@ 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,5 +958,5 @@ 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) Index: /trunk/LMDZ.GENERIC/libf/phystd/phyetat0_mod.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/phyetat0_mod.F90 (revision 3334) +++ /trunk/LMDZ.GENERIC/libf/phystd/phyetat0_mod.F90 (revision 3335) @@ -7,6 +7,6 @@ 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,4 +18,5 @@ 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,4 +49,5 @@ 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,5 +131,5 @@ 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,4 +145,25 @@ write(*,*) "phyetat0: Bare ground albedo 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 " + 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 range:", & + minval(albedo), maxval(albedo) ! ZMEA Index: /trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90 (revision 3334) +++ /trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90 (revision 3335) @@ -134,5 +134,5 @@ 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,4 +141,5 @@ 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,4 +160,5 @@ 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,4 +191,7 @@ 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) Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 (revision 3334) +++ /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 (revision 3335) @@ -39,4 +39,5 @@ 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,6 +593,6 @@ 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,9 +627,6 @@ ! 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,5 +2370,5 @@ 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) Index: /trunk/LMDZ.GENERIC/libf/phystd/surfini.F =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/surfini.F (revision 3334) +++ /trunk/LMDZ.GENERIC/libf/phystd/surfini.F (revision 3335) @@ -1,9 +1,14 @@ + 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,5 +29,5 @@ 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,7 +50,4 @@ 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,17 +57,7 @@ - ! 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,3 +66,5 @@ - END + END SUBROUTINE surfini + + END MODULE surfini_mod