Changeset 5791 for LMDZ6/branches/contrails/libf/phylmd/ocean_slab_mod.f90

Timestamp:

Jul 28, 2025, 7:23:15 PM (7 days ago)

Author:

aborella

Message:

Merge with trunk r5789

Location:

LMDZ6/branches/contrails

Files:

• . (modified) (1 prop)
• libf/phylmd/ocean_slab_mod.f90 (modified) (17 diffs)

LMDZ6/branches/contrails/libf/phylmd/ocean_slab_mod.f90

@@ -50 +50 @@
   !$OMP THREADPRIVATE(fsic)
   ! temperature of the sea ice
-  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
-  !$OMP THREADPRIVATE(tice)
+!GG
+  !REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice_slab
+  !!$OMP THREADPRIVATE(tice_slab)
+  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice_slab
+  !$OMP THREADPRIVATE(tice_slab)
+!GG
   ! sea ice thickness, in kg/m2
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
@@ -238 +242 @@
         ENDIF
         bilg_cum(:) = 0.0   
-        ALLOCATE(tice(klon), stat = error)
+!GG        ALLOCATE(tice(klon), stat = error)
+        ALLOCATE(tice_slab(klon), stat = error)
         IF (error /= 0) THEN
            abort_message='Pb allocation slab_tice'
@@ -633 +638 @@
                       e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
                       ! new ice
-                      tice(ki)=t_freeze
+                !GG      tice(ki)=t_freeze
+                      tice_slab(ki)=t_freeze
                       fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
                       IF (fsic(ki).GT.ice_frac_min) THEN
@@ -650 +656 @@
                       ! quantity of new ice formed over open ocean
                       e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
-                               /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
+                               /(ice_lat+ice_cap/2.*(t_freeze-tice_slab(ki)))
+              !GG                 /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
                       ! new ice height and fraction
                       h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
@@ -759 +766 @@
             cal(i)=2.*RCPD/(snow(i)*ice_cap)
             ! snow conductive flux
-            f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
+            f_cond=sno_cond*(tice_slab(ki)-tsurf_in(i))/snow(i)
+       !GG     f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
             ! all shortwave flux absorbed
             f_swpen=0.
             ! bottom flux (ice conduction)
-            slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
+            slab_bilg(ki)=ice_cond*(tice_slab(ki)-t_freeze)/seaice(ki)
+       !GG     slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
             ! update ice temperature
-            tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
+       !GG     tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
+            tice_slab(ki)=tice_slab(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
                      *(slab_bilg(ki)+f_cond)*dtime
        ELSE ! bare ice
@@ -771 +781 @@
             cal(i)=2.*RCPD/(seaice(ki)*ice_cap)
             ! conductive flux
-            f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
+       !GG     f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
+            f_cond=ice_cond*(t_freeze-tice_slab(ki))/seaice(ki)
             ! penetrative shortwave flux...
             f_swpen=swnet(i)*pen_frac*exp(-pen_ext*seaice(ki)/ice_den)
@@ -789 +800 @@
         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
     DO i=1,knon 
-        IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
+    !GG    IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
+        IF (snow(i).LT.snow_min) tice_slab(knindex(i))=tsurf_new(i)
     END DO
 
@@ -819 +831 @@
             ! energy available for melting snow (in kg of melted snow /m2)
             e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
-               /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
+               /(ice_lat+ice_cap/2.*(t_melt-tice_slab(ki))),0.0),snow(i))
+        !GG       /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
             ! remove snow
             IF (snow(i).GT.e_melt) THEN
@@ -827 +840 @@
                 ! add remaining heat flux to ice
                 e_melt=e_melt-snow(i)
-                tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
-                tsurf_new(i)=tice(ki)
+                tice_slab(ki)=tice_slab(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
+         !GG       tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
+                tsurf_new(i)=tice_slab(ki)
+         !GG       tsurf_new(i)=tice(ki)
             END IF
         END IF
 ! melt ice from above if Tice>0
-        IF (tice(ki).GT.t_melt) THEN
+      !GG  IF (tice(ki).GT.t_melt) THEN
+        IF (tice_slab(ki).GT.t_melt) THEN
             ! quantity of ice melted (kg/m2)
-            e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. & 
+      !GG      e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. & 
+            e_melt=MAX(seaice(ki)*(tice_slab(ki)-t_melt)*ice_cap/2. & 
              /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
             ! melt from above, height only
@@ -850 +867 @@
             seaice(ki)=seaice(ki)-dhsic
             ! surface temperature at melting point
-            tice(ki)=t_melt
+        !GG    tice(ki)=t_melt
+            tice_slab(ki)=t_melt
             tsurf_new(i)=t_melt
         END IF
@@ -860 +878 @@
             seaice(ki)=seaice(ki)+dhsic
             snow(i)=snow(i)-dhsic*sno_den/ice_den
-            ! available energy (freeze sea water + bring to tice)
+            ! available energy (freeze sea water + bring to tice_slab)
             e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ &
-                   ice_cap/2.*(t_freeze-tice(ki)))
+                   ice_cap/2.*(t_freeze-tice_slab(ki)))
+       !GG            ice_cap/2.*(t_freeze-tice(ki)))
             ! update ice temperature
-            tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
+       !GG     tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
+            tice_slab(ki)=tice_slab(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
         END IF
     END DO
@@ -882 +902 @@
         ! ice albedo (varies with ice tkickness and temp)
         alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1)
-        IF (tice(ki).GT.t_freeze-0.01) THEN
+     !GG   IF (tice(ki).GT.t_freeze-0.01) THEN
+        IF (tice_slab(ki).GT.t_freeze-0.01) THEN
             alb_ice=MIN(alb_ice,alb_ice_wet)
         ELSE
@@ -888 +909 @@
         END IF
         ! pond albedo
-        alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
+        alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice_slab(ki)-t_melt,-2.0)))
+     !GG   alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
         ! pond fraction
-        frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
+        frac_pond=0.2*(2.0+MIN(0.0,MAX(tice_slab(ki)-t_melt,-2.0)))
+     !GG   frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
         ! snow fraction
         frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
@@ -917 +940 @@
                ! quantity of new ice
                 e_melt=(t_freeze-tslab(ki,1))/cyang & 
-                       /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
+                       /(ice_lat+ice_cap/2.*(t_freeze-tice_slab(ki)))
+       !GG                /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
                ! first increase height to h_thin
                dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki)))
@@ -935 +959 @@
                ! quantity of melted ice
                e_melt=(tslab(ki,1)-t_freeze)/cyang & 
-                       /(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
+                       /(ice_lat+ice_cap/2.*(tice_slab(ki)-t_freeze))
+          !GG             /(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
                ! first decrease height to h_thick
                dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
@@ -960 +985 @@
     WHERE (fsic.LT.ice_frac_min)
         tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
-        tice=t_melt
+        tice_slab=t_melt
+   !GG     tice=t_melt
         fsic=0.
         seaice=0.
@@ -976 +1002 @@
     IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
     IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
-    IF (ALLOCATED(tice)) DEALLOCATE(tice)
+    IF (ALLOCATED(tice_slab)) DEALLOCATE(tice_slab)
     IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
     IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)