Changeset 12065 for NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/ICE/icethd.F90

Timestamp:

2019-12-05T12:06:36+01:00 (4 years ago)

Author:

smueller

Message:

Synchronizing with ^{/NEMO/trunk@12055 (ticket #2194)}

File:

• NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/ICE/icethd.F90 (modified) (9 diffs)

NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/ICE/icethd.F90

@@ -95 +95 @@
       IF( ln_timing    )   CALL timing_start('icethd')                                                             ! timing
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
 
       IF( kt == nit000 .AND. lwp ) THEN
@@ -102 +103 @@
       ENDIF
       
-      CALL ice_var_glo2eqv
-
       !---------------------------------------------!
       ! computation of friction velocity at T points
@@ -162 +161 @@
             qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
 
-            ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting 
-            IF( zqld > 0._wp ) THEN
+            ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting 
+            ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget
+            IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN
                fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
                qlead(ji,jj) = 0._wp
@@ -178 +178 @@
       ! In case we bypass open-water ice formation
       IF( .NOT. ln_icedO )  qlead(:,:) = 0._wp
-      ! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere
+      ! In case we bypass growing/melting from top and bottom
       IF( .NOT. ln_icedH ) THEN
-         qt_atm_oi  (:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
          qsb_ice_bot(:,:) = 0._wp
          fhld       (:,:) = 0._wp
@@ -221 +220 @@
             dh_i_sub  (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp
             dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp
-            !
-            IF( ln_icedH ) THEN                                     ! --- growing/melting --- !
-                              CALL ice_thd_zdf                             ! Ice/Snow Temperature profile
-                              CALL ice_thd_dh                              ! Ice/Snow thickness   
-                              CALL ice_thd_pnd                             ! Melt ponds formation
-                              CALL ice_thd_ent( e_i_1d(1:npti,:) )         ! Ice enthalpy remapping
+            !                                      
+                              CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- !
+            !
+            IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- !
+                              CALL ice_thd_dh                           ! Ice-Snow thickness   
+                              CALL ice_thd_pnd                          ! Melt ponds formation
+                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping
             ENDIF
-            !
                               CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !    
             !
-                              CALL ice_thd_temp                     ! --- temperature update --- !
+                              CALL ice_thd_temp                     ! --- Temperature update --- !
             !
             IF( ln_icedH .AND. ln_virtual_itd ) &
-               &              CALL ice_thd_mono                     ! --- extra lateral melting if virtual_itd --- !
-            !
-            IF( ln_icedA )    CALL ice_thd_da                       ! --- lateral melting --- !
+               &              CALL ice_thd_mono                     ! --- Extra lateral melting if virtual_itd --- !
+            !
+            IF( ln_icedA )    CALL ice_thd_da                       ! --- Lateral melting --- !
             !
                               CALL ice_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
             !                                                       ! --- & Move to 2D arrays --- !
-            !
          ENDIF
          !
       END DO
-      ! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting)
-      oa_i(:,:,:) = o_i(:,:,:) * a_i(:,:,:)
-
+      !
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
-      !
-                           CALL ice_var_zapsmall           ! --- remove very small ice concentration (<1e-10) --- !
-      !                                                    !     & make sure at_i=SUM(a_i) & ato_i=1 where at_i=0
+      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
       !                   
-      IF( jpl > 1      )   CALL ice_itd_rem( kt )          ! --- Transport ice between thickness categories --- !
-      !
-      IF( ln_icedO     )   CALL ice_thd_do                 ! --- frazil ice growing in leads --- !
+      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- !
+      !
+      IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- !
       !
       ! controls
@@ -418 +412 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m )
          CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m )
-
+         !
+         ! to update ice age
+         CALL tab_2d_1d( npti, nptidx(1:npti), o_i_1d (1:npti), o_i (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
+         !
          ! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
          DO jk = 1, nlay_i
@@ -443 +441 @@
          sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti)
          v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
+         oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti)
          
          CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
@@ -516 +515 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) )
          !
       END SELECT
@@ -541 +541 @@
       REWIND( numnam_ice_ref )              ! Namelist namthd in reference namelist : Ice thermodynamics
       READ  ( numnam_ice_ref, namthd, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd in reference namelist' )
       REWIND( numnam_ice_cfg )              ! Namelist namthd in configuration namelist : Ice thermodynamics
       READ  ( numnam_ice_cfg, namthd, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd in configuration namelist' )
       IF(lwm) WRITE( numoni, namthd )
       !