Index: trunk/LMDZ.GENERIC/changelog.txt =================================================================== --- trunk/LMDZ.GENERIC/changelog.txt (revision 3342) +++ trunk/LMDZ.GENERIC/changelog.txt (revision 3352) @@ -1945,2 +1945,10 @@ Adapts the Thermal Plume Model for generic tracers. Adds variables such as vteta in writediagfi section. + +== 31/05/2024 == EM +Updates for the Slab ocean: +- ensure that surface heat capacity "capcal" at first call matches what it + was at last call of previous run. Also ensure that "ocean slab" variables + put in restartfi.nc are correct. +- adapt "ocean_slab_get_vars" routine to also return sea ice temperature +- add some missing OpenMP THREADPRIVATE on saved variables in ocean_slab_mod. Index: trunk/LMDZ.GENERIC/libf/phystd/ocean_slab_mod.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/ocean_slab_mod.F90 (revision 3342) +++ trunk/LMDZ.GENERIC/libf/phystd/ocean_slab_mod.F90 (revision 3352) @@ -106,4 +106,5 @@ REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice [in W/(m.K) or (W.kg)/(K.m4)] REAL, PRIVATE, SAVE :: sno_cond ! conductivity of snow [in W/(m.K) or (W.kg)/(K.m4)] +!$OMP THREADPRIVATE(sno_cond) REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice [in J/(kg.K)] REAL, PARAMETER :: sea_cap=3994. ! specific heat capacity, seawater [in J/(kg.K)] @@ -116,4 +117,5 @@ REAL, PARAMETER :: ice_frac_min=0.001 REAL, PRIVATE, SAVE :: ice_frac_max ! Max ice fraction (leads) +!$OMP THREADPRIVATE(ice_frac_max) REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness [in kg/m2] REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness [in kg/m2] @@ -121,4 +123,5 @@ ! ice_thin is also height of new ice REAL, PRIVATE, SAVE :: h_ice_thick ! thin ice thickness +!$OMP THREADPRIVATE(h_ice_thick) ! above ice_thick, priority is melt height / grow lateral REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice [in kg/m2] @@ -1049,6 +1052,7 @@ ! ------ End Albedo ---------- - !tests remaining ice fraction - WHERE ((fsic.LT.ice_frac_min).OR.(seaice.LT.h_ice_min)) + !tests remaining ice fraction (on ocean grid points) + WHERE ((zmasq==0).and. & + ((fsic.LT.ice_frac_min).OR.(seaice.LT.h_ice_min))) tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang tice=t_melt @@ -1085,6 +1089,6 @@ !**************************************************************************************** ! - SUBROUTINE ocean_slab_get_vars(ngrid,tslab_loc, seaice_loc, flux_g_loc, & - dt_hdiff_loc,dt_ekman_loc) + SUBROUTINE ocean_slab_get_vars(ngrid, tslab_loc, tice_loc, seaice_loc, & + flux_g_loc, dt_hdiff_loc,dt_ekman_loc) ! "Get some variables from module ocean_slab_mod" @@ -1092,4 +1096,5 @@ INTEGER, INTENT(IN) :: ngrid REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: tslab_loc + REAL, INTENT(OUT) :: tice_loc(ngrid) REAL, DIMENSION(ngrid), INTENT(OUT) :: seaice_loc REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_g_loc @@ -1100,10 +1105,10 @@ ! Set the output variables - tslab_loc(:,:) = 0. + tslab_loc(:,:) = tslab(:,:) + tice_loc(:)=tice(:) + seaice_loc(:) = seaice(:) + flux_g_loc(:) = slab_bilg(:) dt_hdiff_loc(:,:)=0. dt_ekman_loc(:,:)=0. - tslab_loc(:,:) = tslab(:,:) - seaice_loc(:) = seaice(:) - flux_g_loc(:) = slab_bilg(:) !! dt_hdiff_loc(:,i) = dt_hdiff(:,i)*slabh(i)*1000.*4228. !Convert en W/m2 Index: trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 (revision 3342) +++ trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 (revision 3352) @@ -666,5 +666,8 @@ !! CALL init_masquv(ngrid,zmasq) - + ! ensure that "slab" variables match "physiq" ones + call ocean_slab_get_vars(ngrid, tslab, tsea_ice, sea_ice, flux_g, & + dt_hdiff, dt_ekman) + endif ! end of 'ok_slab_ocean'. @@ -678,15 +681,23 @@ if (ok_slab_ocean) then - do ig=1,ngrid - if (nint(rnat(ig)).eq.2) then - capcal(ig)=capcalocean - if (pctsrf_sic(ig).gt.0.5) then - capcal(ig)=capcalseaice - if (qsurf(ig,igcm_h2o_ice).gt.0.) then - capcal(ig)=capcalsno - endif - endif - endif - enddo + ! mimic what is done during a regular time step + do ig=1,ngrid + fluxgrdocean(ig)=fluxgrd(ig) + if (nint(rnat(ig)).eq.0) then + capcal(ig)=capcalocean + fluxgrd(ig)=0. + ! Dividing by cell area to have flux in W/m2 + fluxgrdocean(ig)=flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))/cell_area(ig) + do iq=1,nsoilmx + tsoil(ig,iq)=tsurf(ig) + enddo + if (pctsrf_sic(ig).gt.0.5) then + capcal(ig)=capcalseaice + if (qsurf(ig,igcm_h2o_ice).gt.0.) then + capcal(ig)=capcalsno + endif + endif ! of if (pctsrf_sic(ig).gt.0.5) + endif ! of if (nint(rnat(ig)).eq.0) + enddo ! of do ig=1,ngrid endif ! end of 'ok_slab_ocean'. @@ -1979,5 +1990,5 @@ call ocean_slab_frac(pctsrf_sic, zmasq) - call ocean_slab_get_vars(ngrid, tslab, sea_ice, flux_g, & + call ocean_slab_get_vars(ngrid, tslab, tsea_ice, sea_ice, flux_g, & dt_hdiff, dt_ekman) @@ -2374,4 +2385,8 @@ if (ngrid.ne.1) then write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin + + ! fetch "ocean variables" to ensure they are stored + call ocean_slab_get_vars(ngrid, tslab, tsea_ice, sea_ice, flux_g, & + dt_hdiff, dt_ekman) call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &