Changeset 3352 for trunk/LMDZ.GENERIC/libf

Timestamp:

May 31, 2024, 7:16:42 PM (5 months ago)

Author:

emillour

Message:

Generic PCM:
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.

EM

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• ocean_slab_mod.F90 (modified) (7 diffs)
• physiq_mod.F90 (modified) (4 diffs)

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

@@ -106 +106 @@
     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 +117 @@
     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 +123 @@
     ! 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 +1052 @@
 ! ------ 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 +1089 @@
 !****************************************************************************************
 !
-  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 +1096 @@
     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 +1105 @@
 
 ! 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

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

@@ -666 +666 @@
 
 !!           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 +681 @@
 
             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 +1990 @@
             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 +2385 @@
          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, &