Changeset 15349

Timestamp:

2021-10-08T15:16:52+02:00 (20 months ago)

Author:

cetlod

Message:

dev_PISCO : update with trunk@15348

Location:

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO

Files:

• cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_top_cfg (modified) (3 diffs)
• cfgs/SHARED/field_def_nemo-oce.xml (modified) (1 diff)
• cfgs/SHARED/grid_def_nemo.xml (modified) (1 diff)
• src/ICE/ice.F90 (modified) (3 diffs)
• src/ICE/icealb.F90 (modified) (4 diffs)
• src/ICE/icecor.F90 (modified) (3 diffs)
• src/ICE/icedyn_rdgrft.F90 (modified) (9 diffs)
• src/ICE/icedyn_rhg_evp.F90 (modified) (12 diffs)
• src/ICE/icedyn_rhg_vp.F90 (modified) (15 diffs)
• src/ICE/icethd.F90 (modified) (6 diffs)
• src/ICE/icethd_do.F90 (modified) (9 diffs)
• src/ICE/icethd_pnd.F90 (modified) (2 diffs)
• src/NST/agrif_oce_sponge.F90 (modified) (2 diffs)
• src/NST/agrif_oce_update.F90 (modified) (11 diffs)
• src/NST/agrif_top_update.F90 (modified) (1 diff)
• src/NST/agrif_user.F90 (modified) (7 diffs)
• src/OCE/BDY/bdydyn2d.F90 (modified) (1 diff)
• src/OCE/BDY/bdyini.F90 (modified) (17 diffs)
• src/OCE/DIA/dia25h.F90 (modified) (10 diffs)
• src/OCE/DIA/diahth.F90 (modified) (2 diffs)
• src/OCE/DIA/diawri.F90 (modified) (11 diffs)
• src/OCE/DOM/dom_oce.F90 (modified) (1 diff)
• src/OCE/DOM/domain.F90 (modified) (3 diffs)
• src/OCE/DOM/dommsk.F90 (modified) (1 diff)
• src/OCE/DYN/dynldf_lap_blp_lf.F90 (modified) (1 diff)
• src/OCE/FLO/flodom.F90 (modified) (2 diffs)
• src/OCE/IOM/prtctl.F90 (modified) (2 diffs)
• src/OCE/IOM/restart.F90 (modified) (3 diffs)
• src/OCE/ISF/isfcav.F90 (modified) (1 diff)
• src/OCE/LBC/lbc_lnk_neicoll_generic.h90 (modified) (1 diff)
• src/OCE/LBC/lbc_lnk_pt2pt_generic.h90 (modified) (1 diff)
• src/OCE/LBC/lbc_nfd_generic.h90 (modified) (15 diffs)
• src/OCE/LBC/lbc_nfd_nogather_generic.h90 (deleted)
• src/OCE/LBC/lbcnfd.F90 (modified) (5 diffs)
• src/OCE/LBC/lib_mpp.F90 (modified) (4 diffs)
• src/OCE/LBC/mpp_nfd_generic.h90 (modified) (7 diffs)
• src/OCE/LBC/mppini.F90 (modified) (14 diffs)
• src/OCE/SBC/sbcmod.F90 (modified) (2 diffs)
• src/OCE/SBC/sbcrnf.F90 (modified) (1 diff)
• src/OCE/SBC/sbcssm.F90 (modified) (2 diffs)
• src/OCE/SBC/sbcwave.F90 (modified) (1 diff)
• src/OCE/TRA/eosbn2.F90 (modified) (1 diff)
• src/OCE/TRA/traadv_mus.F90 (modified) (2 diffs)
• src/OCE/USR/usrdef_sbc.F90 (modified) (3 diffs)
• src/OCE/ZDF/zdfevd.F90 (modified) (1 diff)
• src/OCE/ZDF/zdfgls.F90 (modified) (1 diff)
• src/OCE/ZDF/zdfmxl.F90 (modified) (5 diffs)
• src/OCE/ZDF/zdfphy.F90 (modified) (2 diffs)
• src/OCE/ZDF/zdfric.F90 (modified) (1 diff)
• src/OCE/ZDF/zdfsh2.F90 (modified) (1 diff)
• src/OCE/lib_fortran.F90 (modified) (11 diffs)
• src/OCE/lib_fortran_generic.h90 (modified) (2 diffs)
• src/OCE/nemogcm.F90 (modified) (1 diff)
• src/OFF/nemogcm.F90 (modified) (1 diff)
• src/SAO/nemogcm.F90 (modified) (1 diff)
• src/SAS/nemogcm.F90 (modified) (1 diff)
• src/SWE/nemogcm.F90 (modified) (5 diffs)
• src/SWE/stpmlf.F90 (modified) (1 diff)
• src/SWE/stprk3.F90 (modified) (3 diffs)
• src/TOP/PISCES/P4Z/p4zlys.F90 (modified) (2 diffs)
• src/TOP/PISCES/P4Z/p4zsed.F90 (modified) (17 diffs)
• src/TOP/PISCES/P4Z/p4zsms.F90 (modified) (2 diffs)
• src/TOP/PISCES/SED/sedchem.F90 (modified) (21 diffs)
• tests/BENCH/MY_SRC/usrdef_hgr.F90 (modified) (4 diffs)
• tests/BENCH/MY_SRC/usrdef_zgr.F90 (modified) (1 diff)
• tests/DOME/EXPREF/1_namelist_cfg (modified) (1 diff)
• tests/DOME/MY_SRC/usrdef_nam.F90 (modified) (1 diff)
• tests/ISOMIP+/MY_SRC/eosbn2.F90 (modified) (4 diffs)
• tests/OVERFLOW/EXPREF/AGRIF/1_namelist_cfg (modified) (1 diff)
• tests/SWG/MY_SRC/usrdef_sbc.F90 (modified) (1 diff)
• tests/VORTEX/MY_SRC/usrdef_nam.F90 (modified) (1 diff)

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_top_cfg

@@ -65 +65 @@
    sn_trcdta(14) = 'data_FER_nomask.nc',        -1         ,  'Fer'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , ''
    sn_trcdta(23) = 'data_NO3_nomask.nc',        -1         ,  'NO3'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , ''
-   rn_trfac(1)   =   1.028e-06   !  multiplicative factor
-   rn_trfac(2)   =   1.028e-06   !  -      -      -     -
-   rn_trfac(3)   =  44.6e-06     !  -      -      -     -
-   rn_trfac(5)   = 117.0e-06     !  -      -      -     -
-   rn_trfac(7)   =   1.0e-06     !  -      -      -     -
-   rn_trfac(10)  =   1.0e-06     !  -      -      -     -
-   rn_trfac(14)  =   1.0e-06     !  -      -      -     -
-   rn_trfac(23)  =   7.3125e-06  !  -      -      -     -
+   rn_trfac(1)   =   1.028e-06  !  multiplicative factor
+   rn_trfac(2)   =   1.028e-06  !  -      -      -     -
+   rn_trfac(3)   =  44.6e-06  !  -      -      -     -
+   rn_trfac(5)   = 122.0e-06  !  -      -      -     -
+   rn_trfac(7)   =   1.0e-06  !  -      -      -     -
+   rn_trfac(10)  =   1.0e-06  !  -      -      -     -
+   rn_trfac(14)  =   1.0e-06  !  -      -      -     -
+   rn_trfac(23)  =   7.6e-06  !  -      -      -     -
 /
 !-----------------------------------------------------------------------
@@ -118 +118 @@
    rn_trsfac(14) = 6.2667860e-04    !  (  0.035 / 55.85 )
    rn_trsfac(23) = 5.2232143e-01    !  ( From kgN m-2 s-1 to molC l-1 ====> zfact = 7.3125/14 )
-   rn_sbc_time   =  1.              !  Time scaling factor for SBC and CBC data (seconds in a day)
+   rn_sbc_time   = 1.               !  Time scaling factor for SBC and CBC data (seconds in a day)
    !
    sn_trccbc(1)  = 'river.orca'      ,    120            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
@@ -127 +127 @@
    sn_trccbc(14) = 'river.orca'      ,    120            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
    sn_trccbc(23) = 'river.orca'      ,    120            , 'riverdin'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
-   rn_trcfac(1)  = 8.333e+01      !  ( data in Mg/m2/yr : 1e3/12/ryyss)
-   rn_trcfac(2)  = 8.333e+01      !  ( 1e3 /12 )
-   rn_trcfac(5)  = 3.774193e+04   !  ( 1e3 / 31. * 117 )
+   rn_trcfac(1)  = 8.333333e+01      !  ( data in Mg/m2/yr : 1e3/12/ryyss)
+   rn_trcfac(2)  = 8.333333e+01      !  ( 1e3 /12 )
+   rn_trcfac(5)  = 3.774193e+03   !  ( 1e3 / 31. * 117 )
    rn_trcfac(7)  = 3.558719e+01   !  ( 1e3 / 28.1 )
    rn_trcfac(10) = 8.333333e+01   !  ( 1e3 / 12

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/cfgs/SHARED/field_def_nemo-oce.xml

@@ -1243 +1243 @@
 
   <!-- 25h diagnostic output -->
-  <field_group id="25h_grid_T" grid_ref="grid_T_3D" operation="instant">
+  <field_group id="25h_grid_T" grid_ref="grid_T_3D_inner" operation="instant">
     <field id="temper25h"         name="potential temperature 25h mean"    unit="degC" />
     <field id="tempis25h"         name="insitu temperature 25h mean"    unit="degC" />
     <field id="salin25h"          name="salinity 25h mean"                 unit="psu"  />
-    <field id="ssh25h"            name="sea surface height 25h mean"  grid_ref="grid_T_2D"      unit="m"    />
-  </field_group>
-
-  <field_group id="25h_grid_U" grid_ref="grid_U_3D" operation="instant" >
+    <field id="ssh25h"            name="sea surface height 25h mean"  grid_ref="grid_T_2D_inner"      unit="m"    />
+  </field_group>
+
+  <field_group id="25h_grid_U" grid_ref="grid_U_3D_inner" operation="instant" >
     <field id="vozocrtx25h"         name="i current 25h mean"    unit="m/s"   />
   </field_group>
 
-  <field_group id="25h_grid_V" grid_ref="grid_V_3D" operation="instant">
+  <field_group id="25h_grid_V" grid_ref="grid_V_3D_inner" operation="instant">
     <field id="vomecrty25h"         name="j current 25h mean"    unit="m/s"    />
   </field_group>
 
-  <field_group id="25h_grid_W" grid_ref="grid_W_3D" operation="instant">
+  <field_group id="25h_grid_W" grid_ref="grid_W_3D_inner" operation="instant">
     <field id="vovecrtz25h"         name="k current 25h mean"                 unit="m/s"      />
     <field id="avt25h"              name="vertical diffusivity25h mean"       unit="m2/s" />

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/cfgs/SHARED/grid_def_nemo.xml

@@ -87 +87 @@
   </grid>
   <grid id="grid_W_3D_inner" >
-    <domain domain_ref="grid_W" />
-    <axis axis_ref="depthw_inner" />
+    <domain domain_ref="grid_W_inner" />
+    <axis axis_ref="depthw" />
   </grid>
   <!--  -->

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/ice.F90

@@ -147 +147 @@
    !
    !                                     !!** ice-ridging/rafting namelist (namdyn_rdgrft) **
+   LOGICAL,  PUBLIC ::   ln_str_H79       !: ice strength parameterization (Hibler79) (may be used in rheology)
    REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength (also used for landfast param)
+   REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength, Hibler JPO79 (may be used in rheology)
    !
    !                                     !!** ice-rheology namelist (namdyn_rhg) **
@@ -251 +253 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_oce,v_oce     !: surface ocean velocity used in ice dynamics
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ht_i_new        !: ice collection thickness accreted in leads
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fraz_frac       !: fraction of frazil ice accreted at the ice bottom
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   strength        !: ice strength
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   stress1_i, stress2_i, stress12_i   !: 1st, 2nd & diagonal stress tensor element
@@ -453 +456 @@
 
       ii = 1
-      ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) , ht_i_new  (jpi,jpj) , strength(jpi,jpj) ,  &
-         &      stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,                      &
-         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) ,                      &
-         &      aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv  (jpi,jpj) , STAT=ierr(ii) )
+      ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) , ht_i_new (jpi,jpj) , fraz_frac (jpi,jpj) ,  &
+         &      strength (jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,  &
+         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i  (jpi,jpj) ,                        &
+         &      aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv (jpi,jpj) , STAT=ierr(ii) )
 
       ii = ii + 1

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icealb.F90

@@ -30 +30 @@
    PUBLIC   ice_alb        ! called in icesbc.F90 and iceupdate.F90
 
-   REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
+   REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066_wp   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
    !
    !                             !!* albedo namelist (namalb)
@@ -111 +111 @@
       REAL(wp) ::   zalb_snw, zafrac_snw      ! snow-covered sea ice albedo & relative snow fraction
       REAL(wp) ::   zalb_cs, zalb_os          ! albedo of ice under clear/overcast sky
+      !! clem
+      REAL(wp), PARAMETER ::   zhi_albcst = 1.5_wp ! pivotal thickness (should be in the namelist)
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('icealb')
       !
-      z1_href_pnd = 1. / 0.05
-      z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
-      z1_c2 = 1. / 0.05
-      z1_c3 = 1. / 0.02
-      z1_c4 = 1. / 0.03
+      z1_href_pnd = 1._wp / 0.05_wp
+      z1_c1 = 1._wp / ( LOG(zhi_albcst) - LOG(0.05_wp) ) 
+      z1_c2 = 1._wp / 0.05_wp
+      z1_c3 = 1._wp / 0.02_wp
+      z1_c4 = 1._wp / 0.03_wp
       !
       CALL ice_var_snwfra( ph_snw, za_s_fra )   ! calculate ice fraction covered by snow
@@ -148 +150 @@
             ENDIF
             !                       !--- Bare ice albedo (for hi < 150cm)
-            IF( 0.05 < ph_ice(ji,jj,jl) .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN      ! 5cm < hi < 150cm
-               zalb_ice = zalb_ice    + ( 0.18 - zalb_ice   ) * z1_c1 * ( LOG(1.5) - LOG(ph_ice(ji,jj,jl)) )
-            ELSEIF( ph_ice(ji,jj,jl) <= 0.05 ) THEN                               ! 0cm < hi < 5cm
-               zalb_ice = rn_alb_oce  + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice(ji,jj,jl)
+            IF( 0.05 < ph_ice(ji,jj,jl) .AND. ph_ice(ji,jj,jl) <= zhi_albcst ) THEN      ! 5cm < hi < 150cm
+               zalb_ice = zalb_ice    + ( 0.18_wp - zalb_ice   ) * z1_c1 * ( LOG(zhi_albcst) - LOG(ph_ice(ji,jj,jl)) )
+            ELSEIF( ph_ice(ji,jj,jl) <= 0.05_wp ) THEN                               ! 0cm < hi < 5cm
+               zalb_ice = rn_alb_oce  + ( 0.18_wp - rn_alb_oce ) * z1_c2 * ph_ice(ji,jj,jl)
             ENDIF
             !
@@ -166 +168 @@
             zalb_os = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1)
             !
-            zalb_cs = zalb_os - ( - 0.1010 * zalb_os * zalb_os  &
-               &                  + 0.1933 * zalb_os - 0.0148 ) * tmask(ji,jj,1)
+            zalb_cs = zalb_os - ( - 0.1010_wp * zalb_os * zalb_os  &
+               &                  + 0.1933_wp * zalb_os - 0.0148_wp ) * tmask(ji,jj,1)
             !
             ! albedo depends on cloud fraction because of non-linear spectral effects

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icecor.F90

@@ -53 +53 @@
       INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo
       !
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
+      INTEGER  ::   ji, jj, jl       ! dummy loop indices
       REAL(wp) ::   zsal, zzc
       !!----------------------------------------------------------------------
@@ -99 +99 @@
          END DO
       ENDIF
-
+      !
       IF( kn /= 0 ) THEN   ! no zapsmall if kn=0 (for bdy for instance) because we do not want ice-ocean exchanges (wfx,sfx,hfx)
          !                                                              otherwise conservation diags will fail
@@ -105 +105 @@
          CALL ice_var_zapsmall      !  Zap small values                                  !
          !                          !-----------------------------------------------------
-      ENDIF
-      !                             !-----------------------------------------------------
-      IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values !
-         DO_2D( 0, 0, 0, 0 )        !-----------------------------------------------------
-            IF ( at_i(ji,jj) == 0._wp ) THEN    ! what to do if there is no ice
-               IF ( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side
-               IF ( at_i(ji-1,jj) == 0._wp )   u_ice(ji-1,jj) = 0._wp   ! left side
-               IF ( at_i(ji,jj+1) == 0._wp )   v_ice(ji,jj  ) = 0._wp   ! upper side
-               IF ( at_i(ji,jj-1) == 0._wp )   v_ice(ji,jj-1) = 0._wp   ! bottom side
-            ENDIF
-         END_2D
-         CALL lbc_lnk( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp )
       ENDIF
       !

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icedyn_rdgrft.F90

@@ -57 +57 @@
    !
    ! ** namelist (namdyn_rdgrft) **
-   LOGICAL  ::   ln_str_H79       ! ice strength parameterization (Hibler79)
-   REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79
    REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging
    LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975))
@@ -272 +270 @@
 
       ! controls
-      IF( sn_cfctl%l_prtctl )   CALL ice_prt3D   ('icedyn_rdgrft')                                                        ! prints
+      IF( sn_cfctl%l_prtctl )   CALL ice_prt3D('icedyn_rdgrft')                                                           ! prints
       IF( ln_icectl    )   CALL ice_prt     (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ')                             ! prints
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
@@ -520 +518 @@
       !
       INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation
+      LOGICAL , DIMENSION(jpij) ::   ll_shift         ! logical for doing calculation or not
       !!-------------------------------------------------------------------
       !
@@ -540 +539 @@
          DO ji = 1, npti
 
-            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
+            ! set logical to true when ridging
+            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ;   ll_shift(ji) = .TRUE.
+            ELSE                                                                ;   ll_shift(ji) = .FALSE.
+            ENDIF
+            
+            IF( ll_shift(ji) ) THEN   ! only if ice is ridging
 
                IF( a_i_2d(ji,jl1) > epsi10 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
@@ -630 +634 @@
          DO jk = 1, nlay_s
             DO ji = 1, npti
-               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
+               IF( ll_shift(ji) ) THEN
                   ! Compute ridging /rafting fractions
                   afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
@@ -651 +655 @@
          DO jk = 1, nlay_i
             DO ji = 1, npti
-               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
+               IF( ll_shift(ji) ) THEN
                   ! Compute ridging /rafting fractions
                   afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji)
@@ -674 +678 @@
             DO ji = 1, npti
 
-               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
+               IF( ll_shift(ji) ) THEN
 
                   ! Compute the fraction of ridged ice area and volume going to thickness category jl2
@@ -731 +735 @@
             DO jk = 1, nlay_s
                DO ji = 1, npti
-                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
+                  IF( ll_shift(ji) )   &
                      &   ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji)  +  &
                      &                                               esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
@@ -740 +744 @@
             DO jk = 1, nlay_i
                DO ji = 1, npti
-                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
+                  IF( ll_shift(ji) )   &
                      &   ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)
                END DO

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icedyn_rhg_evp.F90

@@ -52 +52 @@
    INTEGER ::   nvarid   ! netcdf variable id
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   fimask   ! mask at F points for the ice
-
+   
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -180 +180 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)
+      !! -- advect fields at the rheology time step for the calculation of strength
+      LOGICAL  ::   ll_advups = .FALSE.
+      REAL(wp) ::   zdt_ups
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   ztmp
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   za_i_ups, zv_i_ups   ! tracers advected upstream
       !!-------------------------------------------------------------------
 
@@ -234 +239 @@
          z1_alph2 = 1._wp / ( zalph2 + 1._wp )
       ELSE
-         zdtevp   = rdt_ice
+         zdtevp   = rDt_ice
          ! zalpha parameters set later on adaptatively
       ENDIF
@@ -280 +285 @@
 
          ! Ocean currents at U-V points
-         v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
-         u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
+         ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
+         v_oceU(ji,jj)   = 0.25_wp * ( (v_oce(ji,jj) + v_oce(ji,jj-1)) + (v_oce(ji+1,jj) + v_oce(ji+1,jj-1)) ) * umask(ji,jj,1)
+         u_oceV(ji,jj)   = 0.25_wp * ( (u_oce(ji,jj) + u_oce(ji-1,jj)) + (u_oce(ji,jj+1) + u_oce(ji-1,jj+1)) ) * vmask(ji,jj,1)
 
          ! m/dt
@@ -382 +388 @@
             zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
 
-         END_2D
-         CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp )
-
-         ! P/delta at T points
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            ! P/delta at T points
             zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
-         END_2D
-
+
+         END_2D
+         CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp, zp_delt, 'T', 1.0_wp )
+
+         !
          DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )   ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
 
             ! divergence at T points (duplication to avoid communications)
-            zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-               &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+            ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
+            zdiv  = ( (e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj))   &
+               &    + (e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1))   &
                &    ) * r1_e1e2t(ji,jj)
 
@@ -442 +448 @@
 
             ! P/delta at F points
-            zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
+            ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
+            zp_delf = 0.25_wp * ( (zp_delt(ji,jj) + zp_delt(ji+1,jj)) + (zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1)) )
 
             ! stress at F points (zkt/=0 if landfast)
@@ -450 +457 @@
 
          ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
+         ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             !                   !--- U points
-            zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
+            zfU(ji,jj) = 0.5_wp * ( (( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
                &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
-               &                    ) * r1_e2u(ji,jj)                                                                      &
+               &                    ) * r1_e2u(ji,jj))                                                                      &
                &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
                &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
@@ -460 +468 @@
             !
             !                !--- V points
-            zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
+            zfV(ji,jj) = 0.5_wp * ( (( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
                &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
-               &                    ) * r1_e1v(ji,jj)                                                                      &
+               &                    ) * r1_e1v(ji,jj))                                                                      &
                &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
                &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
@@ -468 +476 @@
             !
             !                !--- ice currents at U-V point
-            v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1)
-            u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1)
+            v_iceU(ji,jj) = 0.25_wp * ( (v_ice(ji,jj) + v_ice(ji,jj-1)) + (v_ice(ji+1,jj) + v_ice(ji+1,jj-1)) ) * umask(ji,jj,1)
+            u_iceV(ji,jj) = 0.25_wp * ( (u_ice(ji,jj) + u_ice(ji-1,jj)) + (u_ice(ji,jj+1) + u_ice(ji-1,jj+1)) ) * vmask(ji,jj,1)
             !
          END_2D
@@ -682 +690 @@
          IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 )
          !
+         !
+         ! --- change strength according to advected a_i and v_i (upstream for now) --- !
+         IF( ll_advups .AND. ln_str_H79 ) THEN
+            !
+            IF( jter == 1 ) THEN                               ! init
+               ALLOCATE( za_i_ups(jpi,jpj,jpl), zv_i_ups(jpi,jpj,jpl) )
+               ALLOCATE( ztmp(jpi,jpj) )
+               zdt_ups = rDt_ice / REAL( nn_nevp )
+               za_i_ups(:,:,:) = a_i(:,:,:)
+               zv_i_ups(:,:,:) = v_i(:,:,:)
+            ELSE
+               CALL lbc_lnk( 'icedyn_rhg_evp', za_i_ups, 'T', 1.0_wp, zv_i_ups, 'T', 1.0_wp )               
+            ENDIF
+            !
+            CALL rhg_upstream( zdt_ups, u_ice, v_ice, za_i_ups )   ! upstream advection: a_i
+            CALL rhg_upstream( zdt_ups, u_ice, v_ice, zv_i_ups )   ! upstream advection: v_i
+            !
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )    ! strength
+               strength(ji,jj) = rn_pstar * SUM( zv_i_ups(ji,jj,:) ) * EXP( -rn_crhg * ( 1._wp - SUM( za_i_ups(ji,jj,:) ) ) )
+            END_2D
+            IF( nn_hls == 1 )  CALL lbc_lnk( 'icedyn_rhg_evp', strength, 'T', 1.0_wp )
+            !
+            DO_2D( 0, 0, 0, 0 )                                ! strength smoothing
+               IF( SUM( za_i_ups(ji,jj,:) ) > 0._wp ) THEN
+                  ztmp(ji,jj) = ( 4._wp * strength(ji,jj) + strength(ji-1,jj  ) + strength(ji+1,jj  ) &
+                     &                                    + strength(ji  ,jj-1) + strength(ji  ,jj+1) &
+                     &          ) / ( 4._wp + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
+               ELSE
+                  ztmp(ji,jj) = 0._wp
+               ENDIF
+            END_2D
+            DO_2D( 0, 0, 0, 0 )
+               strength(ji,jj) = ztmp(ji,jj)
+            END_2D
+            !
+            IF( jter == nn_nevp ) THEN
+               DEALLOCATE( za_i_ups, zv_i_ups )
+               DEALLOCATE( ztmp )
+            ENDIF
+         ENDIF
          !                                                ! ==================== !
       END DO                                              !  end loop over jter  !
       !                                                   ! ==================== !
       IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta )
+      !
+      IF( ll_advups .AND. ln_str_H79 )   CALL lbc_lnk( 'icedyn_rhg_evp', strength, 'T', 1.0_wp )
       !
       !------------------------------------------------------------------------------!
@@ -760 +810 @@
       IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
       IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
+      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
 
       ! --- Stress tensor invariants (SIMIP diags) --- !
@@ -1018 +1069 @@
    END SUBROUTINE rhg_evp_rst
 
+   SUBROUTINE rhg_upstream( pdt, pu, pv, pt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE rhg_upstream  ***
+      !!
+      !! **  Purpose :   compute the upstream fluxes and upstream guess of tracer
+      !!----------------------------------------------------------------------
+      REAL(wp)                   , INTENT(in   ) ::   pdt              ! tracer time-step
+      REAL(wp), DIMENSION(:,:  ) , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
+      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   pt               ! tracer fields
+      !
+      INTEGER  ::   ji, jj, jl    ! dummy loop indices
+      REAL(wp) ::   ztra          ! local scalar
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups   ! upstream fluxes
+      !!----------------------------------------------------------------------
+      DO jl = 1, jpl
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
+            zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji  ,jj  ,jl) + &
+               &             MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji+1,jj  ,jl)
+            zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji  ,jj  ,jl) + &
+               &             MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji  ,jj+1,jl)
+         END_2D
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) )
+            !
+            pt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
+         END_2D
+      END DO
+      !
+   END SUBROUTINE rhg_upstream
 
 #else

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icedyn_rhg_vp.F90

@@ -218 +218 @@
       
       !!----------------------------------------------------------------------------------------------------------------------
-      ! DEBUG put all forcing terms to zero
-         ! air-ice drag
-         utau_ice(:,:) = 0._wp
-         vtau_ice(:,:) = 0._wp
-         ! coriolis
-         ff_t(:,:) = 0._wp
-         ! ice-ocean drag
-         rn_cio = 0._wp
-         ! ssh 
-         ! done line 330 !!! dont forget to act there
-      ! END DEBUG
+!!$      ! DEBUG put all forcing terms to zero
+!!$         ! air-ice drag
+!!$         utau_ice(:,:) = 0._wp
+!!$         vtau_ice(:,:) = 0._wp
+!!$         ! coriolis
+!!$         ff_t(:,:) = 0._wp
+!!$         ! ice-ocean drag
+!!$         rn_cio = 0._wp
+!!$         ! ssh 
+!!$         ! done line 330 !!! dont forget to act there
+!!$      ! END DEBUG
 
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)'
@@ -240 +240 @@
       
       ! for diagnostics and convergence tests
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
       END_2D
       IF( nn_rhg_chkcvg > 0 ) THEN
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
          END_2D
@@ -329 +329 @@
       !    non-embedded sea ice: use ocean surface for slope calculation
       zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
-      zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!!
-
-      zmt(:,:) = rhos * vt_s(:,:) + rhoi * vt_i(:,:)       ! Snow and ice mass at T-point
-      zmf(:,:) = zmt(:,:) * ff_t(:,:)                      ! Coriolis factor at T points (m*f)
-      
-      DO jj = 2, jpj - 1
-         DO ji = 2, jpi - 1
-
-            ! Ice fraction at U-V points
-            za_iU(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
-            za_iV(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
-
-            ! Snow and ice mass at U-V points
-            zm1             = zmt(ji,jj)
-            zm2             = zmt(ji+1,jj)
-            zm3             = zmt(ji,jj+1)
-            
-            zmassU          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
-            zmassV          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
-                          
-            ! Mass per unit area divided by time step
-            zmassU_t(ji,jj) = zmassU * r1_Dt_ice
-            zmassV_t(ji,jj) = zmassV * r1_Dt_ice
-
-            ! Acceleration term contribution to RHS (depends on velocity at previous time step)            
-            zmU_t(ji,jj)    = zmassU_t(ji,jj) * u_ice(ji,jj)
-            zmV_t(ji,jj)    = zmassV_t(ji,jj) * v_ice(ji,jj)
-            
-            ! Ocean currents at U-V points
-            v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
-            u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
-            
-            ! Wind stress
-            ztaux_ai(ji,jj) = za_iU(ji,jj) * utau_ice(ji,jj)
-            ztauy_ai(ji,jj) = za_iV(ji,jj) * vtau_ice(ji,jj)
-
-            ! Force due to sea surface tilt(- m*g*GRAD(ssh))
-            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
-            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
-
-            ! Mask for ice presence (1) or absence (0)
-            zmsk00x(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
-            zmsk00y(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
-
-            ! Mask for lots of ice (1) or little ice (0)
-            IF ( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
-            ELSE                                                      ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
-            IF ( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
-            ELSE                                                      ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF              
+!!$      zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!!
+
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         zm1          = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) )  ! Ice/snow mass at U-V points
+         zmf  (ji,jj) = zm1 * ff_t(ji,jj)                            ! Coriolis at T points (m*f)
+      END_2D
+      
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+
+         ! Ice fraction at U-V points
+         za_iU(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+         za_iV(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+         
+         ! Ice/snow mass at U-V points
+         zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
+         zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
+         zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )         
+         zmassU          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+         zmassV          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+         
+         ! Mass per unit area divided by time step
+         zmassU_t(ji,jj) = zmassU * r1_Dt_ice
+         zmassV_t(ji,jj) = zmassV * r1_Dt_ice
+         
+         ! Acceleration term contribution to RHS (depends on velocity at previous time step)            
+         zmU_t(ji,jj)    = zmassU_t(ji,jj) * u_ice(ji,jj)
+         zmV_t(ji,jj)    = zmassV_t(ji,jj) * v_ice(ji,jj)
+         
+         ! Ocean currents at U-V points
+         v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
+         u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
+         
+         ! Wind stress
+         ztaux_ai(ji,jj) = za_iU(ji,jj) * utau_ice(ji,jj)
+         ztauy_ai(ji,jj) = za_iV(ji,jj) * vtau_ice(ji,jj)
+         
+         ! Force due to sea surface tilt(- m*g*GRAD(ssh))
+         zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
+         zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
+         
+         ! Mask for ice presence (1) or absence (0)
+         zmsk00x(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
+         zmsk00y(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
+         
+         ! Mask for lots of ice (1) or little ice (0)
+         IF ( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
+         ELSE                                                      ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
+         IF ( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
+         ELSE                                                      ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF              
 
 ! MV TEST DEBUG
-            IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji+1,jj)  <= zmmin )     .AND.  &
-               & ( at_i(ji,jj)  <= zamin .OR. at_i(ji+1,jj) <= zamin ) )    THEN   ;   zmsk01x(ji,jj) = 0._wp
-            ELSE                                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
-
-            IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji,jj+1)  <= zmmin )     .AND.  &
-               & ( at_i(ji,jj)  <= zamin .OR. at_i(ji,jj+1) <= zamin ) )    THEN   ;   zmsk01y(ji,jj) = 0._wp
-            ELSE                                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF              
+!!$            IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji+1,jj)  <= zmmin )     .AND.  &
+!!$               & ( at_i(ji,jj)  <= zamin .OR. at_i(ji+1,jj) <= zamin ) )    THEN   ;   zmsk01x(ji,jj) = 0._wp
+!!$            ELSE                                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
+!!$
+!!$            IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji,jj+1)  <= zmmin )     .AND.  &
+!!$               & ( at_i(ji,jj)  <= zamin .OR. at_i(ji,jj+1) <= zamin ) )    THEN   ;   zmsk01y(ji,jj) = 0._wp
+!!$            ELSE                                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF              
 ! END MV TEST DEBUG
 
-         END DO
-      END DO   
-
-      CALL iom_put( 'zmsk00x'    , zmsk00x  )   ! MV DEBUG
-      CALL iom_put( 'zmsk00y'    , zmsk00y  )   ! MV DEBUG
-      CALL iom_put( 'zmsk01x'    , zmsk01x  )   ! MV DEBUG
-      CALL iom_put( 'zmsk01y'    , zmsk01y  )   ! MV DEBUG
-      CALL iom_put( 'ztaux_ai'   , ztaux_ai )   ! MV DEBUG
-      CALL iom_put( 'ztauy_ai'   , ztauy_ai )   ! MV DEBUG
-      CALL iom_put( 'zspgU'      , zspgU    )   ! MV DEBUG
-      CALL iom_put( 'zspgV'      , zspgV    )   ! MV DEBUG
-            
+      END_2D
+
+!!$      CALL iom_put( 'zmsk00x'    , zmsk00x  )   ! MV DEBUG
+!!$      CALL iom_put( 'zmsk00y'    , zmsk00y  )   ! MV DEBUG
+!!$      CALL iom_put( 'zmsk01x'    , zmsk01x  )   ! MV DEBUG
+!!$      CALL iom_put( 'zmsk01y'    , zmsk01y  )   ! MV DEBUG
+!!$      CALL iom_put( 'ztaux_ai'   , ztaux_ai )   ! MV DEBUG
+!!$      CALL iom_put( 'ztauy_ai'   , ztauy_ai )   ! MV DEBUG
+!!$      CALL iom_put( 'zspgU'      , zspgU    )   ! MV DEBUG
+!!$      CALL iom_put( 'zspgV'      , zspgV    )   ! MV DEBUG
+!!$            
       !------------------------------------------------------------------------------!
       !
@@ -437 +436 @@
 
          ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
-         DO jj = 1, jpj - 1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
-            DO ji = 1, jpi - 1
-
-               ! shear at F points
-               zds(ji,jj) = ( ( zu_c(ji,jj+1) * r1_e1u(ji,jj+1) - zu_c(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
-                  &         + ( zv_c(ji+1,jj) * r1_e2v(ji+1,jj) - zv_c(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-                  &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
-
-            END DO
-         END DO
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
+
+            ! shear at F points
+            zds(ji,jj) = ( ( zu_c(ji,jj+1) * r1_e1u(ji,jj+1) - zu_c(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+               &         + ( zv_c(ji+1,jj) * r1_e2v(ji+1,jj) - zv_c(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+               &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
+
+         END_2D
 
          CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan
-         CALL iom_put( 'zds'          , zds      )   ! MV DEBUG
+!!$         CALL iom_put( 'zds'          , zds      )   ! MV DEBUG
 
          IF( lwp )   WRITE(numout,*) ' outer loop  1a i_out : ', i_out
@@ -496 +493 @@
          CALL lbc_lnk( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. )
 
-         CALL iom_put( 'zzt'        , zzt      )   ! MV DEBUG
-         CALL iom_put( 'zet'        , zet      )   ! MV DEBUG
-         CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG
+!!$         CALL iom_put( 'zzt'        , zzt      )   ! MV DEBUG
+!!$         CALL iom_put( 'zet'        , zet      )   ! MV DEBUG
+!!$         CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG
 
          IF( lwp )   WRITE(numout,*) ' outer loop  1b i_out : ', i_out
@@ -515 +512 @@
          
          CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. )
-         CALL iom_put( 'zef'          , zef            ) ! MV DEBUG
+!!$         CALL iom_put( 'zef'          , zef            ) ! MV DEBUG
          IF( lwp )   WRITE(numout,*) ' outer loop  1c i_out : ', i_out
 
@@ -558 +555 @@
          CALL lbc_lnk( 'icedyn_rhg_vp', zCorU,  'U', -1., zCorV, 'V', -1. )
 
-         CALL iom_put( 'zCwU'          , zCwU           ) ! MV DEBUG
-         CALL iom_put( 'zCwV'          , zCwV           ) ! MV DEBUG
-         CALL iom_put( 'zCorU'         , zCorU          ) ! MV DEBUG
-         CALL iom_put( 'zCorV'         , zCorV          ) ! MV DEBUG
+!!$         CALL iom_put( 'zCwU'          , zCwU           ) ! MV DEBUG
+!!$         CALL iom_put( 'zCwV'          , zCwV           ) ! MV DEBUG
+!!$         CALL iom_put( 'zCorU'         , zCorU          ) ! MV DEBUG
+!!$         CALL iom_put( 'zCorV'         , zCorV          ) ! MV DEBUG
 
          IF( lwp )   WRITE(numout,*) ' outer loop  1f i_out : ', i_out
@@ -591 +588 @@
          END DO
 
-         CALL iom_put( 'zs1_rhsu'      , zs1_rhsu       ) ! MV DEBUG
-         CALL iom_put( 'zs2_rhsu'      , zs2_rhsu       ) ! MV DEBUG
-         CALL iom_put( 'zs1_rhsv'      , zs1_rhsv       ) ! MV DEBUG
-         CALL iom_put( 'zs2_rhsv'      , zs2_rhsv       ) ! MV DEBUG
+!!$         CALL iom_put( 'zs1_rhsu'      , zs1_rhsu       ) ! MV DEBUG
+!!$         CALL iom_put( 'zs2_rhsu'      , zs2_rhsu       ) ! MV DEBUG
+!!$         CALL iom_put( 'zs1_rhsv'      , zs1_rhsv       ) ! MV DEBUG
+!!$         CALL iom_put( 'zs2_rhsv'      , zs2_rhsv       ) ! MV DEBUG
          
          ! a priori, no lbclnk, because rhsu is only used in the inner domain
@@ -619 +616 @@
          CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. )
 
-         CALL iom_put( 'zs12_rhsu'     , zs12_rhsu      ) ! MV DEBUG
-         CALL iom_put( 'zs12_rhsv'     , zs12_rhsv      ) ! MV DEBUG
+!!$         CALL iom_put( 'zs12_rhsu'     , zs12_rhsu      ) ! MV DEBUG
+!!$         CALL iom_put( 'zs12_rhsv'     , zs12_rhsv      ) ! MV DEBUG
 
          ! a priori, no lbclnk, because rhsu are only used in the inner domain
@@ -644 +641 @@
          END DO
 
-         CALL iom_put( 'zf_rhsu'       , zf_rhsu        ) ! MV DEBUG
-         CALL iom_put( 'zf_rhsv'       , zf_rhsv        ) ! MV DEBUG
+!!$         CALL iom_put( 'zf_rhsu'       , zf_rhsu        ) ! MV DEBUG
+!!$         CALL iom_put( 'zf_rhsv'       , zf_rhsv        ) ! MV DEBUG
          
          !---------------------------
@@ -666 +663 @@
          CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V',  -1.)
 
-         CALL iom_put( 'zmU_t'         , zmU_t          ) ! MV DEBUG
-         CALL iom_put( 'zmV_t'         , zmV_t          ) ! MV DEBUG
-         CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu  ) ! MV DEBUG
-         CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv  ) ! MV DEBUG
-         CALL iom_put( 'zrhsu'         , zrhsu          ) ! MV DEBUG
-         CALL iom_put( 'zrhsv'         , zrhsv          ) ! MV DEBUG
-         
+!!$         CALL iom_put( 'zmU_t'         , zmU_t          ) ! MV DEBUG
+!!$         CALL iom_put( 'zmV_t'         , zmV_t          ) ! MV DEBUG
+!!$         CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu  ) ! MV DEBUG
+!!$         CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv  ) ! MV DEBUG
+!!$         CALL iom_put( 'zrhsu'         , zrhsu          ) ! MV DEBUG
+!!$         CALL iom_put( 'zrhsv'         , zrhsv          ) ! MV DEBUG
+!!$         
          ! inner domain calculations -> no lbclnk
 
@@ -773 +770 @@
          CALL lbc_lnk( 'icedyn_rhg_vp', zEU  , 'U', 1., zEV  , 'V',  1. )
                
-         CALL iom_put( 'zAU'           , zAU            ) ! MV DEBUG
-         CALL iom_put( 'zBU'           , zBU            ) ! MV DEBUG
-         CALL iom_put( 'zCU'           , zCU            ) ! MV DEBUG
-         CALL iom_put( 'zDU'           , zDU            ) ! MV DEBUG
-         CALL iom_put( 'zEU'           , zEU            ) ! MV DEBUG
-         CALL iom_put( 'zAV'           , zAV            ) ! MV DEBUG
-         CALL iom_put( 'zBV'           , zBV            ) ! MV DEBUG
-         CALL iom_put( 'zCV'           , zCV            ) ! MV DEBUG
-         CALL iom_put( 'zDV'           , zDV            ) ! MV DEBUG
-         CALL iom_put( 'zEV'           , zEV            ) ! MV DEBUG
+!!$         CALL iom_put( 'zAU'           , zAU            ) ! MV DEBUG
+!!$         CALL iom_put( 'zBU'           , zBU            ) ! MV DEBUG
+!!$         CALL iom_put( 'zCU'           , zCU            ) ! MV DEBUG
+!!$         CALL iom_put( 'zDU'           , zDU            ) ! MV DEBUG
+!!$         CALL iom_put( 'zEU'           , zEU            ) ! MV DEBUG
+!!$         CALL iom_put( 'zAV'           , zAV            ) ! MV DEBUG
+!!$         CALL iom_put( 'zBV'           , zBV            ) ! MV DEBUG
+!!$         CALL iom_put( 'zCV'           , zCV            ) ! MV DEBUG
+!!$         CALL iom_put( 'zDV'           , zDV            ) ! MV DEBUG
+!!$         CALL iom_put( 'zEV'           , zEV            ) ! MV DEBUG
 
       !------------------------------------------------------------------------------!
@@ -1103 +1100 @@
       CALL lbc_lnk( 'icedyn_rhg_vp', zCU_prime  , 'U',  1., zCV_prime  , 'V',  1. )
 
-      CALL iom_put( 'zFU'           , zFU            ) ! MV DEBUG
-      CALL iom_put( 'zBU_prime'     , zBU_prime      ) ! MV DEBUG
-      CALL iom_put( 'zCU_prime'     , zCU_prime      ) ! MV DEBUG
-      CALL iom_put( 'zFU_prime'     , zFU_prime      ) ! MV DEBUG
-
-      CALL iom_put( 'zFV'           , zFV            ) ! MV DEBUG
-      CALL iom_put( 'zBV_prime'     , zBV_prime      ) ! MV DEBUG
-      CALL iom_put( 'zCV_prime'     , zCV_prime      ) ! MV DEBUG
-      CALL iom_put( 'zFV_prime'     , zFV_prime      ) ! MV DEBUG
+!!$      CALL iom_put( 'zFU'           , zFU            ) ! MV DEBUG
+!!$      CALL iom_put( 'zBU_prime'     , zBU_prime      ) ! MV DEBUG
+!!$      CALL iom_put( 'zCU_prime'     , zCU_prime      ) ! MV DEBUG
+!!$      CALL iom_put( 'zFU_prime'     , zFU_prime      ) ! MV DEBUG
+!!$
+!!$      CALL iom_put( 'zFV'           , zFV            ) ! MV DEBUG
+!!$      CALL iom_put( 'zBV_prime'     , zBV_prime      ) ! MV DEBUG
+!!$      CALL iom_put( 'zCV_prime'     , zCV_prime      ) ! MV DEBUG
+!!$      CALL iom_put( 'zFV_prime'     , zFV_prime      ) ! MV DEBUG
 
       CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
 
-      IF ( lwp ) WRITE(numout,*) ' We are about to output uice_dbg '
-      IF( iom_use('uice_dbg' ) )   CALL iom_put( 'uice_dbg'   , u_ice    )                              ! ice velocity u after solver
-      IF( iom_use('vice_dbg' ) )   CALL iom_put( 'vice_dbg'   , v_ice    )                              ! ice velocity v after solver
+!!$      IF ( lwp ) WRITE(numout,*) ' We are about to output uice_dbg '
+!!$      IF( iom_use('uice_dbg' ) )   CALL iom_put( 'uice_dbg'   , u_ice    )                              ! ice velocity u after solver
+!!$      IF( iom_use('vice_dbg' ) )   CALL iom_put( 'vice_dbg'   , v_ice    )                              ! ice velocity v after solver
             
       !------------------------------------------------------------------------------!
@@ -1134 +1131 @@
       ENDIF
 
-      ! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead
-      DO jj = 2, jpj - 1
-         DO ji = 2, jpi - 1   
-
-             u_ice(ji,jj) =   zmsk00x(ji,jj)                               & 
-      &         * (           zmsk01x(ji,jj)   * u_oce(ji,jj) * 0.01_wp    &
-                  + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    )
-
-             v_ice(ji,jj) =   zmsk00y(ji,jj)                               & 
-      &         * (           zmsk01y(ji,jj)   * v_oce(ji,jj) * 0.01_wp    &
-                  + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    )
-
-         END DO
-      END DO
-
-      CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
-
-      IF ( lwp ) WRITE(numout,*) ' Velocity replaced '
+!!$      ! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead
+!!$      DO jj = 2, jpj - 1
+!!$         DO ji = 2, jpi - 1   
+!!$
+!!$             u_ice(ji,jj) =   zmsk00x(ji,jj)                               & 
+!!$      &         * (           zmsk01x(ji,jj)   * u_oce(ji,jj) * 0.01_wp    &
+!!$                  + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    )
+!!$
+!!$             v_ice(ji,jj) =   zmsk00y(ji,jj)                               & 
+!!$      &         * (           zmsk01y(ji,jj)   * v_oce(ji,jj) * 0.01_wp    &
+!!$                  + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    )
+!!$
+!!$         END DO
+!!$      END DO
+!!$
+!!$      CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
+!!$
+!!$      IF ( lwp ) WRITE(numout,*) ' Velocity replaced '
 
       ! END DEBUG
@@ -1365 +1362 @@
                zfac             =   zp_deltastar_t(ji,jj)
                zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) )
-               zsig1            = 0._wp !!! FUCKING DEBUG TEST !!!
+!!$               zsig1            = 0._wp !!! FUCKING DEBUG TEST !!!
                zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
                zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icethd.F90

@@ -20 +20 @@
    USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl
    USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, &
-      &                 qml_ice, qcn_ice, qtr_ice_top
+      &                 qml_ice, qcn_ice, qtr_ice_top, utau_ice, vtau_ice
    USE ice1D          ! sea-ice: thermodynamics variables
    USE icethd_zdf     ! sea-ice: vertical heat diffusion
@@ -97 +97 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric, zvel   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
       !
+      ! for collection thickness
+      INTEGER  ::   iter             !   -       -
+      REAL(wp) ::   zvfrx, zvgx, ztaux, zf                     !   -      -
+      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp   !   -      -
+      REAL(wp), PARAMETER ::   zcai = 1.4e-3_wp               ! ice-air drag (clem: should be dependent on coupling/forcing used)
+      REAL(wp), PARAMETER ::   zhicrit = 0.04_wp                                             ! frazil ice thickness
+      REAL(wp), PARAMETER ::   zsqcd   = 1.0_wp / SQRT( 1.3_wp * zcai )                      ! 1/SQRT(airdensity*drag)
+      REAL(wp), PARAMETER ::   zgamafr = 0.03_wp
       !!-------------------------------------------------------------------
       ! controls
@@ -128 +136 @@
                &                         ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) )
          END_2D
-      ELSE      !  if no ice dynamics => transmit directly the atmospheric stress to the ocean
+      ELSE      !  if no ice dynamics => transfer directly the atmospheric stress to the ocean
          DO_2D( 0, 0, 0, 0 )
             zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp *  &
@@ -219 +227 @@
       ENDIF
 
+
+      !---------------------------------------------------------!
+      ! Collection thickness of ice formed in leads and polynyas
+      !---------------------------------------------------------!    
+      ! ht_i_new is the thickness of new ice formed in open water
+      ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T)
+      ! Frazil ice forms in open water, is transported by wind, accumulates at the edge of the consolidated ice edge
+      ! where it forms aggregates of a specific thickness called collection thickness.
+      !
+      fraz_frac(:,:) = 0._wp
+      !
+      ! Default new ice thickness
+      WHERE( qlead(:,:) < 0._wp ) ! cooling
+         ht_i_new(:,:) = rn_hinew
+      ELSEWHERE
+         ht_i_new(:,:) = 0._wp
+      END WHERE
+
+      IF( ln_frazil ) THEN
+         ztwogp  = 2._wp * rho0 / ( grav * 0.3_wp * ( rho0 - rhoi ) )  ! reduced grav
+         !
+         DO_2D( 0, 0, 0, 0 )
+            IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling
+               ! -- Wind stress -- !
+               ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   &
+                  &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp
+               ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   &
+                  &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp
+               ! Square root of wind stress
+               ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
+
+               ! -- Frazil ice velocity -- !
+               rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
+               zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
+               zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
+
+               ! -- Pack ice velocity -- !
+               zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
+               zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
+
+               ! -- Relative frazil/pack ice velocity -- !
+               rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
+               zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
+                  &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15_wp * 0.15_wp ) * rswitch
+
+               !!clem
+               fraz_frac(ji,jj) = rswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz
+               !!clem
+               
+               ! -- new ice thickness (iterative loop) -- !
+               ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1_wp )    &
+                  &                      / ( ( zhicrit + 0.1_wp ) * ( zhicrit + 0.1_wp ) -  zhicrit * zhicrit ) * ztwogp * zvrel2
+               iter = 1
+               DO WHILE ( iter < 20 ) 
+                  zf  = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) -   &
+                     &    ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2
+                  zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0_wp * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
+
+                  ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 )
+                  iter = iter + 1
+               END DO
+               !
+               ! bound ht_i_new (though I don't see why it should be necessary)
+               ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) )
+               !
+            ELSE
+               ht_i_new(ji,jj) = 0._wp
+            ENDIF
+            !
+         END_2D
+         ! 
+         CALL lbc_lnk( 'icethd', fraz_frac, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp  )
+
+      ENDIF
+
       !-------------------------------------------------------------------------------------------!
       ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories
@@ -268 +351 @@
       !
       IF ( ln_pnd .AND. ln_icedH ) &
-         &                    CALL ice_thd_pnd                      ! --- Melt ponds
+         &                    CALL ice_thd_pnd                      ! --- Melt ponds --- !
       !
       IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- !
@@ -276 +359 @@
                               CALL ice_cor( kt , 2 )                ! --- Corrections --- !
       !
-      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice              ! ice natural aging incrementation
+      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice              ! --- Ice natural aging incrementation
+      !
+      DO_2D( 0, 0, 0, 0 )                                           ! --- Ice velocity corrections
+         IF( at_i(ji,jj) == 0._wp ) THEN   ! if ice has melted
+            IF( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side
+            IF( at_i(ji-1,jj) == 0._wp )   u_ice(ji-1,jj) = 0._wp   ! left side
+            IF( at_i(ji,jj+1) == 0._wp )   v_ice(ji,jj  ) = 0._wp   ! upper side
+            IF( at_i(ji,jj-1) == 0._wp )   v_ice(ji,jj-1) = 0._wp   ! bottom side
+         ENDIF
+      END_2D
+      CALL lbc_lnk( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp )
       !
       ! convergence tests

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icethd_do.F90

@@ -17 +17 @@
    USE phycst         ! physical constants
    USE sbc_oce , ONLY : sss_m
-   USE sbc_ice , ONLY : utau_ice, vtau_ice
    USE ice1D          ! sea-ice: thermodynamics variables
    USE ice            ! sea-ice: variables
@@ -38 +37 @@
 
    !                          !!** namelist (namthd_do) **
-   REAL(wp) ::   rn_hinew      ! thickness for new ice formation (m)
-   LOGICAL  ::   ln_frazil     ! use of frazil ice collection as function of wind (T) or not (F)
-   REAL(wp) ::   rn_maxfraz    ! maximum portion of frazil ice collecting at the ice bottom
-   REAL(wp) ::   rn_vfraz      ! threshold drift speed for collection of bottom frazil ice
-   REAL(wp) ::   rn_Cfraz      ! squeezing coefficient for collection of bottom frazil ice
+   REAL(wp), PUBLIC ::   rn_hinew      ! thickness for new ice formation (m)
+   LOGICAL , PUBLIC ::   ln_frazil     ! use of frazil ice collection as function of wind (T) or not (F)
+   REAL(wp), PUBLIC ::   rn_maxfraz    ! maximum portion of frazil ice collecting at the ice bottom
+   REAL(wp), PUBLIC ::   rn_vfraz      ! threshold drift speed for collection of bottom frazil ice
+   REAL(wp), PUBLIC ::   rn_Cfraz      ! squeezing coefficient for collection of bottom frazil ice
 
    !! * Substitutions
@@ -78 +77 @@
       !!------------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      INTEGER  ::   iter             !   -       -
-      REAL(wp) ::   ztmelts, zfrazb, zweight, zde                               ! local scalars
-      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
-      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
-      !
+      !
+      REAL(wp) ::   ztmelts
+      REAL(wp) ::   zdE
       REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
       REAL(wp) ::   zEi          ! sea ice specific enthalpy (J/kg)
@@ -102 +99 @@
       REAL(wp), DIMENSION(jpij) ::   zda_res     ! residual area in case of excessive heat budget
       REAL(wp), DIMENSION(jpij) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
-      REAL(wp), DIMENSION(jpij) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
+      REAL(wp), DIMENSION(jpij) ::   zfraz_frac_1d ! relative ice / frazil velocity (1D vector)
       !
       REAL(wp), DIMENSION(jpij,jpl) ::   zv_b    ! old volume of ice in category jl
@@ -109 +106 @@
       REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d !: 1-D version of e_i
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::   zvrel    ! relative ice / frazil velocity
-      !
-      REAL(wp) :: zcai = 1.4e-3_wp               ! ice-air drag (clem: should be dependent on coupling/forcing used)
       !!-----------------------------------------------------------------------!
 
@@ -119 +113 @@
       at_i(:,:) = SUM( a_i, dim=3 )
       !------------------------------------------------------------------------------!
-      ! 1) Collection thickness of ice formed in leads and polynyas
-      !------------------------------------------------------------------------------!    
-      ! ht_i_new is the thickness of new ice formed in open water
-      ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T)
-      ! Frazil ice forms in open water, is transported by wind
-      ! accumulates at the edge of the consolidated ice edge
-      ! where it forms aggregates of a specific thickness called
-      ! collection thickness.
-
-      zvrel(:,:) = 0._wp
-
-      ! Default new ice thickness
-      WHERE( qlead(:,:) < 0._wp ) ! cooling
-         ht_i_new(:,:) = rn_hinew
-      ELSEWHERE
-         ht_i_new(:,:) = 0._wp
-      END WHERE
-
-      IF( ln_frazil ) THEN
-         !
-         ht_i_new(:,:) = 0._wp
-         !
-         ! Physical constants
-         zhicrit = 0.04                                          ! frazil ice thickness
-         ztwogp  = 2. * rho0 / ( grav * 0.3 * ( rho0 - rhoi ) )  ! reduced grav
-         zsqcd   = 1.0 / SQRT( 1.3 * zcai )                      ! 1/SQRT(airdensity*drag)
-         zgamafr = 0.03
-         !
-         DO_2D( 0, 0, 0, 0 )
-            IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling
-               ! -- Wind stress -- !
-               ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   &
-                  &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp
-               ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   &
-                  &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp
-               ! Square root of wind stress
-               ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
-
-               ! -- Frazil ice velocity -- !
-               rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
-               zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
-               zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
-
-               ! -- Pack ice velocity -- !
-               zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
-               zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
-
-               ! -- Relative frazil/pack ice velocity -- !
-               rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
-               zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
-                  &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch
-               zvrel(ji,jj) = SQRT( zvrel2 )
-
-               ! -- new ice thickness (iterative loop) -- !
-               ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    &
-                  &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2
-
-               iter = 1
-               DO WHILE ( iter < 20 ) 
-                  zf  = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) -   &
-                     &    ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2
-                  zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0 * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
-
-                  ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 )
-                  iter = iter + 1
-               END DO
-               !
-               ! bound ht_i_new (though I don't see why it should be necessary)
-               ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) )
-               !
-            ENDIF
-            !
-         END_2D
-         ! 
-         CALL lbc_lnk( 'icethd_do', zvrel, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp  )
-
-      ENDIF
-
-      !------------------------------------------------------------------------------!
-      ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice
+      ! 1) Compute thickness, salinity, enthalpy, age, area and volume of new ice
       !------------------------------------------------------------------------------!
       ! it occurs if cooling
@@ -223 +138 @@
             END DO
          END DO
-         CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d  (1:npti) , qlead      )
-         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo       )
-         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw    )
-         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw    )
-         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new   )
-         CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d  (1:npti) , zvrel      )
+         CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d     (1:npti) , qlead      )
+         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d      (1:npti) , t_bo       )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d   (1:npti) , sfx_opw    )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d   (1:npti) , wfx_opw    )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice    (1:npti) , ht_i_new   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zfraz_frac_1d(1:npti) , fraz_frac  )
 
          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd    )
@@ -300 +215 @@
          END DO
          
-         zv_frazb(1:npti) = 0._wp
-         IF( ln_frazil ) THEN
-            ! A fraction zfrazb of frazil ice is accreted at the ice bottom
-            DO ji = 1, npti
-               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
-               zfrazb        = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz
-               zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
-               zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
-            END DO
-         END IF
+         ! A fraction fraz_frac of frazil ice is accreted at the ice bottom
+         DO ji = 1, npti
+            rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
+            zv_frazb(ji)  =           zfraz_frac_1d(ji) * rswitch   * zv_newice(ji)
+            zv_newice(ji) = ( 1._wp - zfraz_frac_1d(ji) * rswitch ) * zv_newice(ji)
+         END DO
          
          ! --- Area of new ice --- !
@@ -317 +228 @@
 
          !------------------------------------------------------------------------------!
-         ! 3) Redistribute new ice area and volume into ice categories                  !
+         ! 2) Redistribute new ice area and volume into ice categories                  !
          !------------------------------------------------------------------------------!
 

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icethd_pnd.F90

@@ -662 +662 @@
                      IF ( t_su(ji,jj,jl) > zTp ) THEN
 
-                        zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) )
+                        zdvice = MIN( -dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) )
 
                         IF ( zdvice > epsi10 ) THEN
@@ -775 +775 @@
                ! recalculate equivalent pond variables
                IF ( a_ip(ji,jj,jl) > epsi10) THEN
-                  h_ip(ji,jj,jl)      = v_ip(ji,jj,jl) / a_i(ji,jj,jl)
+                  h_ip(ji,jj,jl)      = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
                   a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar
                   h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/NST/agrif_oce_sponge.F90

@@ -688 +688 @@
             END DO
 
-            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
+            ubdiff(i1:i2,j1:j2,1:jpk) = (uu(i1:i2,j1:j2,1:jpk,Kbb_a) - tabres_child(i1:i2,j1:j2,1:jpk))*umask(i1:i2,j1:j2,1:jpk)
          ELSE
 
-            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
+            ubdiff(i1:i2,j1:j2,1:jpk) = (uu(i1:i2,j1:j2,1:jpk,Kbb_a) - tabres(i1:i2,j1:j2,1:jpk,1))*umask(i1:i2,j1:j2,1:jpk)
   
          ENDIF
@@ -872 +872 @@
             END DO
 
-            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)  
+            vbdiff(i1:i2,j1:j2,1:jpk) = (vv(i1:i2,j1:j2,1:jpk,Kbb_a) - tabres_child(i1:i2,j1:j2,1:jpk))*vmask(i1:i2,j1:j2,1:jpk)  
          ELSE
 
-            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
+            vbdiff(i1:i2,j1:j2,1:jpk) = (vv(i1:i2,j1:j2,1:jpk,Kbb_a) - tabres(i1:i2,j1:j2,1:jpk,1))*vmask(i1:i2,j1:j2,1:jpk)
 
          ENDIF

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/NST/agrif_oce_update.F90

@@ -135 +135 @@
       IF (Agrif_Root()) RETURN
       !
-      Agrif_UseSpecialValueInUpdate = .TRUE.
+      l_vremap = ln_vert_remap
+      Agrif_UseSpecialValueInUpdate = .NOT.l_vremap 
       Agrif_SpecialValueFineGrid = 0._wp
 # if ! defined DECAL_FEEDBACK_2D
@@ -144 +145 @@
       !
       Agrif_UseSpecialValueInUpdate = .FALSE.
+      l_vremap = .FALSE.
       !
 #  if defined VOL_REFLUX
@@ -173 +175 @@
       !       
       Agrif_UseSpecialValueInUpdate = .TRUE.
-      Agrif_SpecialValueFineGrid = 0.
+      Agrif_SpecialValueFineGrid = 0._wp
 
       CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN  )
@@ -199 +201 @@
       !
 #else
-      Agrif_UseSpecialValueInUpdate = .TRUE.
-      Agrif_SpecialValueFineGrid = 0.
+      Agrif_UseSpecialValueInUpdate = .FALSE.
+      Agrif_SpecialValueFineGrid = 0._wp
       ! 
       ! No interface separation here, update vertical grid at T points 
@@ -271 +273 @@
          ! Update total depths:
          ! --------------------
-      hu(:,:,Kmm_a) = 0._wp                        ! Ocean depth at U-points
-      hv(:,:,Kmm_a) = 0._wp                        ! Ocean depth at V-points
+      hu(:,:,Kmm_a) = 0._wp                    ! Ocean depth at U-points
+      hv(:,:,Kmm_a) = 0._wp                    ! Ocean depth at V-points
       DO jk = 1, jpkm1
          hu(:,:,Kmm_a) = hu(:,:,Kmm_a) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk)
@@ -546 +548 @@
 
          ELSE
-            DO jk=1,jpk
+            DO jk=k1,k2
                DO jj=j1,j2
                   DO ji=i1,i2
@@ -657 +659 @@
                   N_out = 0
                   DO jk=1,jpk
-                     IF (vmask(ji,jj,jk) == 0) EXIT
+                     IF (vmask(ji,jj,jk) == 0._wp) EXIT
                      N_out = N_out + 1
                      h_out(N_out) = e3v(ji,jj,jk,Kmm_a)
@@ -693 +695 @@
 
          ELSE
-            DO jk=1,jpk
+            DO jk=k1,k2
                DO jj=j1,j2
                   DO ji=i1,i2
@@ -1097 +1099 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN 
+                  IF( tabres(ji,jj,jk,1) .NE. 0._wp ) THEN 
                      print *,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
                      print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk)
@@ -1226 +1228 @@
                      &                                        ( e3t(ji,jj,jk  ,Kbb_a) - e3t_0(ji,jj,jk  ) )
                      gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a)
-                     gdept(ji,jj,jk,Kbb_a) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a))  &
+                     gdept(ji,jj,jk,Kbb_a) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb_a) + 0.5_wp * e3w(ji,jj,jk,Kbb_a))  &
                          &               + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) +       e3w(ji,jj,jk,Kbb_a)) 
                   END DO
@@ -1259 +1261 @@
                &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t(ji,jj,jk  ,Kmm_a) - e3t_0(ji,jj,jk  ) )
                gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a)
-               gdept(ji,jj,jk,Kmm_a) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a))  &
+               gdept(ji,jj,jk,Kmm_a) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm_a) + 0.5_wp * e3w(ji,jj,jk,Kmm_a))  &
                    &               + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) +       e3w(ji,jj,jk,Kmm_a)) 
                gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/NST/agrif_top_update.F90

@@ -120 +120 @@
                   N_out = 0
                   DO jk=1,jpk ! jpk of parent grid
-                     IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF
+                     IF (tmask(ji,jj,jk) == 0._wp ) EXIT ! TODO: Will not work with ISF
                      N_out = N_out + 1
                      h_out(N_out) = e3t(ji,jj,jk,Kmm_a)

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/NST/agrif_user.F90

@@ -434 +434 @@
       ! 1. Declaration of the type of variable which have to be interpolated
       !---------------------------------------------------------------------
-      ind1 =              nbghostcells
+      ind1 =              nbghostcells - 1 ! Remove one land cell in ghosts 
       ind2 = nn_hls + 1 + nbghostcells_x_w
       ind3 = nn_hls + 1 + nbghostcells_y_s
@@ -634 +634 @@
       !                            2,2 = two ghost lines
       !-------------------------------------------------------------------------------------
-      ind1 =              nbghostcells
+      ind1 =              nbghostcells - 1 ! Remove one land cell in ghosts 
       ind2 = nn_hls + 1 + nbghostcells_x_w
       ind3 = nn_hls + 1 + nbghostcells_y_s
@@ -661 +661 @@
       ! 3. Set location of interpolations
       !----------------------------------
-      CALL Agrif_Set_bc(tra_ice_id,(/0,ind1/))
-      CALL Agrif_Set_bc(u_ice_id  ,(/0,ind1/))
-      CALL Agrif_Set_bc(v_ice_id  ,(/0,ind1/))
-
-      CALL Agrif_Set_bc(tra_iceini_id,(/0,ind1/))
-      CALL Agrif_Set_bc(u_iceini_id  ,(/0,ind1/))
-      CALL Agrif_Set_bc(v_iceini_id  ,(/0,ind1/))
+      CALL Agrif_Set_bc(tra_ice_id,(/0,ind1-1/))
+      CALL Agrif_Set_bc(u_ice_id  ,(/0,ind1-1/))
+      CALL Agrif_Set_bc(v_ice_id  ,(/0,ind1-1/))
+
+      CALL Agrif_Set_bc(tra_iceini_id,(/0,ind1-1/))
+      CALL Agrif_Set_bc(u_iceini_id  ,(/0,ind1-1/))
+      CALL Agrif_Set_bc(v_iceini_id  ,(/0,ind1-1/))
 
       ! 4. Set update type in case 2 ways (child=>parent) (normal & tangent to the grid cell for velocities)
@@ -782 +782 @@
       ! 1. Declaration of the type of variable which have to be interpolated
       !---------------------------------------------------------------------
-      ind1 =              nbghostcells
+      ind1 =              nbghostcells - 1 ! Remove one land cell in ghosts 
       ind2 = nn_hls + 1 + nbghostcells_x_w
       ind3 = nn_hls + 1 + nbghostcells_y_s
@@ -838 +838 @@
       !
       INTEGER  ::   ios                 ! Local integer output status for namelist read
+      INTEGER  ::   imin, imax, jmin, jmax
       NAMELIST/namagrif/ ln_agrif_2way, ln_init_chfrpar, rn_sponge_tra, rn_sponge_dyn, rn_trelax_tra, rn_trelax_dyn, &
                        & ln_spc_dyn, ln_vert_remap, ln_chk_bathy
       !!--------------------------------------------------------------------------------------
       !
-      READ  ( numnam_ref, namagrif, IOSTAT = ios, ERR = 901)
+      IF ( .NOT.Agrif_Root() ) THEN
+         !
+         READ  ( numnam_ref, namagrif, IOSTAT = ios, ERR = 901)
 901 IF( ios /= 0 )   CALL ctl_nam ( ios , 'namagrif in reference namelist' )
-      READ  ( numnam_cfg, namagrif, IOSTAT = ios, ERR = 902 )
+         READ  ( numnam_cfg, namagrif, IOSTAT = ios, ERR = 902 )
 902 IF( ios >  0 )   CALL ctl_nam ( ios , 'namagrif in configuration namelist' )
-      IF(lwm) WRITE ( numond, namagrif )
-      !
-      IF(lwp) THEN                    ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'agrif_nemo_init : AGRIF parameters'
-         WRITE(numout,*) '~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namagrif : set AGRIF parameters'
-         WRITE(numout,*) '      Two way nesting activated ln_agrif_2way         = ', ln_agrif_2way
-         WRITE(numout,*) '      child initial state from parent ln_init_chfrpar = ', ln_init_chfrpar
-         WRITE(numout,*) '      ad. sponge coeft for tracers      rn_sponge_tra = ', rn_sponge_tra
-         WRITE(numout,*) '      ad. sponge coeft for dynamics     rn_sponge_tra = ', rn_sponge_dyn
-         WRITE(numout,*) '      ad. time relaxation for tracers   rn_trelax_tra = ', rn_trelax_tra
-         WRITE(numout,*) '      ad. time relaxation for dynamics  rn_trelax_dyn = ', rn_trelax_dyn
-         WRITE(numout,*) '      use special values for dynamics   ln_spc_dyn    = ', ln_spc_dyn
-         WRITE(numout,*) '      vertical remapping                ln_vert_remap = ', ln_vert_remap
-         WRITE(numout,*) '      check bathymetry                  ln_chk_bathy  = ', ln_chk_bathy
-      ENDIF
-
-! JC => side effects of lines below to be checked:
-      IF (.not.agrif_root()) THEN
-         nbghostcells_x_w = nbghostcells
-         nbghostcells_x_e = nbghostcells
-         nbghostcells_y_s = nbghostcells
-         nbghostcells_y_n = nbghostcells
- 
-
-         lk_west  = .TRUE.
-         lk_east  = .TRUE.
-         lk_south = .TRUE.
-         lk_north = .TRUE.
+         IF(lwm) WRITE ( numond, namagrif )
          !
-         ! Correct number of ghost cells according to periodicity
-         !
-         IF( l_Iperio         ) THEN ; lk_west  = .FALSE. ; lk_east = .FALSE. ; nbghostcells_x_w = 0 ; nbghostcells_x_e = 0 ; ENDIF
-         IF( l_NFold          ) THEN 
-            lk_west   = .FALSE. ; lk_east  = .FALSE. ; nbghostcells_x_w = 0 ; nbghostcells_x_e = 0 
-            lk_north  = .FALSE. ; lk_south = .FALSE. ; nbghostcells_y_s = 0 ; nbghostcells_y_n = 0 
-         ENDIF
-         IF( Agrif_Iy() == 1  ) THEN ; lk_south = .FALSE. ; nbghostcells_y_s = 1 ; ENDIF
-         IF( Agrif_Iy() + nbcellsy/AGRIF_Irhoy() ==  Agrif_Parent(Nj0glo) - 1 ) THEN ; lk_north = .FALSE. ; nbghostcells_y_n = 1 ; ENDIF
-         IF( Agrif_Ix() == 1  ) THEN ; lk_west  = .FALSE. ; nbghostcells_x_w = 1 ; ENDIF
-         IF( Agrif_Ix() + nbcellsx/AGRIF_Irhox() ==  Agrif_Parent(Ni0glo) - 1 ) THEN ; lk_east  = .FALSE. ; nbghostcells_x_e = 1 ; ENDIF
-         !
-         ! Some checks
-         IF( (.NOT.ln_vert_remap).AND.(jpkglo>Agrif_Parent(jpkglo)) )                    CALL ctl_stop( 'STOP',    &
-           &   'agrif_nemo_init: Agrif children must have less or equal number of vertical levels without ln_vert_remap defined' ) 
-         IF( Ni0glo /= nbcellsx + nbghostcells_x_w + nbghostcells_x_e   )   CALL ctl_stop( 'STOP',    &
-           &   'agrif_nemo_init: Agrif children requires jpiglo == nbcellsx + nbghostcells_x_w + nbghostcells_x_e' )
-         IF( Nj0glo /= nbcellsy + nbghostcells_y_s + nbghostcells_y_n )   CALL ctl_stop( 'STOP',    &
-           &   'agrif_nemo_init: Agrif children requires jpjglo == nbcellsy + nbghostcells_y_s + nbghostcells_y_n' )
-         IF( ln_use_jattr )   CALL ctl_stop( 'STOP', 'agrif_nemo_init:Agrif children requires ln_use_jattr = .false. ' )
-      ELSE
-         ! Root grid
-         nbghostcells_x_w = 1 
-         nbghostcells_x_e = 1 
-         nbghostcells_y_s = 1 
-         nbghostcells_y_n = 1 
-         IF ( l_Iperio.OR.l_NFold ) THEN
-           nbghostcells_x_w = 0
-           nbghostcells_x_e = 0
-         ENDIF
-         IF ( l_NFold ) THEN
-           nbghostcells_y_n = 0 ! for completeness
-         ENDIF
-      ENDIF
+         IF(lwp) THEN                    ! control print
+            WRITE(numout,*)
+            WRITE(numout,*) 'agrif_nemo_init : AGRIF parameters'
+            WRITE(numout,*) '~~~~~~~~~~~~~~~'
+            WRITE(numout,*) '   Namelist namagrif : set AGRIF parameters'
+            WRITE(numout,*) '      Two way nesting activated ln_agrif_2way         = ', ln_agrif_2way
+            WRITE(numout,*) '      child initial state from parent ln_init_chfrpar = ', ln_init_chfrpar
+            WRITE(numout,*) '      ad. sponge coeft for tracers      rn_sponge_tra = ', rn_sponge_tra
+            WRITE(numout,*) '      ad. sponge coeft for dynamics     rn_sponge_tra = ', rn_sponge_dyn
+            WRITE(numout,*) '      ad. time relaxation for tracers   rn_trelax_tra = ', rn_trelax_tra
+            WRITE(numout,*) '      ad. time relaxation for dynamics  rn_trelax_dyn = ', rn_trelax_dyn
+            WRITE(numout,*) '      use special values for dynamics   ln_spc_dyn    = ', ln_spc_dyn
+            WRITE(numout,*) '      vertical remapping                ln_vert_remap = ', ln_vert_remap
+            WRITE(numout,*) '      check bathymetry                  ln_chk_bathy  = ', ln_chk_bathy
+         ENDIF
+
+         imin = Agrif_Ix()
+         imax = Agrif_Ix() + nbcellsx/AGRIF_Irhox()
+         jmin = Agrif_Iy()
+         jmax = Agrif_Iy() + nbcellsy/AGRIF_Irhoy()
+         lk_west  = .TRUE. ; lk_east  = .TRUE.
+         lk_north = .TRUE. ; lk_south = .TRUE.
+
+         ! Check zoom position along i:
+         ! ----------------------------
+         IF ( imin >= imax ) THEN
+            CALL ctl_stop( 'STOP', 'AGRIF zoom imin must be < imax' )
+         ENDIF
+
+         IF ( Agrif_Parent(l_Iperio) ) THEN
+            IF ( l_Iperio ) THEN ! Cyclic east-west zoom
+               lk_west = .FALSE. ; lk_east = .FALSE.
+               ! Checks:
+               IF ( imin/=1-Agrif_Parent(nbghostcells_x_w) ) THEN
+                  WRITE(ctmp1, 9000) ' AGRIF zoom is East-West cyclic, imin must = ', &
+                  1 - Agrif_Parent(nbghostcells_x_w)
+                  CALL ctl_stop( 'STOP', ctmp1 )
+               ENDIF
+               IF ( imax/=Agrif_Parent(Ni0glo)+1-Agrif_Parent(nbghostcells_x_w)) THEN
+                  WRITE(ctmp1, 9000) ' AGRIF zoom is East-West cyclic, imax must = ', &
+                  Agrif_Parent(Ni0glo) + 1 - Agrif_Parent(nbghostcells_x_w)
+                  CALL ctl_stop( 'STOP', ctmp1 )
+               ENDIF
+            ELSE
+               IF ( imax>Agrif_Parent(Ni0glo)-Agrif_Parent(nbghostcells_x_w)) THEN
+                  WRITE(ctmp1, 9000) ' AGRIF zoom imax must be <= ', &
+                  Agrif_Parent(Ni0glo) - Agrif_Parent(nbghostcells_x_w)
+                  CALL ctl_stop( 'STOP', ctmp1 )
+               ENDIF
+            ENDIF
+         ELSE
+            IF ( imin<2-Agrif_Parent(nbghostcells_x_w) ) THEN
+               WRITE(ctmp1, 9000) ' AGRIF zoom imin must be >= ', &
+               2 - Agrif_Parent(nbghostcells_x_w)
+               CALL ctl_stop( 'STOP', ctmp1 )
+            ENDIF
+            IF ( imax>Agrif_Parent(Ni0glo)-Agrif_Parent(nbghostcells_x_w)) THEN
+               WRITE(ctmp1, 9000) ' AGRIF zoom imax must be <= ', &
+               Agrif_Parent(Ni0glo) - Agrif_Parent(nbghostcells_x_w)
+               CALL ctl_stop( 'STOP', ctmp1 )
+            ENDIF
+            IF ( imin==2-Agrif_Parent(nbghostcells_x_w) )                    lk_west = .FALSE.
+            IF ( imax==Agrif_Parent(Ni0glo)-Agrif_Parent(nbghostcells_x_w) ) lk_east = .FALSE.  
+         ENDIF
+
+         ! Check zoom position along j:
+         ! ----------------------------
+         IF ( jmin >= jmax ) THEN
+            CALL ctl_stop( 'STOP', 'AGRIF zoom jmin must be < jmax' )
+         ENDIF
+
+         IF ( Agrif_Parent(l_NFold) ) THEN
+            IF ( l_NFold ) THEN ! North-Fold 
+               lk_north = .FALSE.
+               ! Checks:
+               IF ( jmax/=Agrif_Parent(Nj0glo)+1-Agrif_Parent(nbghostcells_y_s)) THEN 
+                  WRITE(ctmp1, 9000) ' AGRIF zoom has a North-Fold, jmax must = ', &
+                  Agrif_Parent(Nj0glo) + 1 - Agrif_Parent(nbghostcells_y_s)
+                  CALL ctl_stop( 'STOP', ctmp1 )
+               ENDIF
+            ENDIF
+         ELSE
+            IF ( jmax>Agrif_Parent(Nj0glo)-Agrif_Parent(nbghostcells_y_s)) THEN 
+               WRITE(ctmp1, 9000) ' AGRIF zoom jmax must be <= ', &
+               Agrif_Parent(Nj0glo) - Agrif_Parent(nbghostcells_y_s)
+               CALL ctl_stop( 'STOP', ctmp1 )
+            ENDIF
+            IF ( jmax==Agrif_Parent(Nj0glo)-Agrif_Parent(nbghostcells_y_s) ) lk_north = .FALSE. 
+         ENDIF
+
+         IF ( jmin<2-Agrif_Parent(nbghostcells_y_s)) THEN 
+            WRITE(ctmp1, 9000) ' AGRIF zoom jmin must be >= ', &
+            2 - Agrif_Parent(nbghostcells_y_s)
+            CALL ctl_stop( 'STOP', ctmp1 )
+         ENDIF
+         IF ( jmin==2-Agrif_Parent(nbghostcells_y_s) ) lk_south = .FALSE. 
+
+      ELSE ! Root grid
+         lk_west  = .FALSE. ; lk_east  = .FALSE.
+         lk_north = .FALSE. ; lk_south = .FALSE.
+      ENDIF
+  
+      ! Set ghost cells including over Parent grid: 
+      nbghostcells_x_w = nbghostcells
+      nbghostcells_x_e = nbghostcells
+      nbghostcells_y_s = nbghostcells
+      nbghostcells_y_n = nbghostcells
+
+      IF (.NOT.lk_west ) nbghostcells_x_w = 1
+      IF (.NOT.lk_east ) nbghostcells_x_e = 1
+      IF (.NOT.lk_south) nbghostcells_y_s = 1
+      IF (.NOT.lk_north) nbghostcells_y_n = 1
+
+      IF ( l_Iperio ) THEN
+         nbghostcells_x_w = 0 ; nbghostcells_x_e = 0
+      ENDIF
+      IF ( l_NFold ) THEN
+         nbghostcells_y_n = 0
+      ENDIF
+      
+      IF ( .NOT.Agrif_Root() ) THEN ! Check expected grid size: 
+         IF( (.NOT.ln_vert_remap).AND.(jpkglo>Agrif_Parent(jpkglo)) )   CALL ctl_stop( 'STOP',    &
+           &   'AGRIF children must have less or equal number of vertical levels without ln_vert_remap defined' ) 
+         IF( Ni0glo /= nbcellsx + nbghostcells_x_w + nbghostcells_x_e ) CALL ctl_stop( 'STOP',    &
+           &   'AGRIF children requires jpiglo == nbcellsx + nbghostcells_x_w + nbghostcells_x_e' )
+         IF( Nj0glo /= nbcellsy + nbghostcells_y_s + nbghostcells_y_n ) CALL ctl_stop( 'STOP',    &
+           &   'AGRIF children requires jpjglo == nbcellsy + nbghostcells_y_s + nbghostcells_y_n' )
+         IF( ln_use_jattr )   CALL ctl_stop( 'STOP', 'AGRIF children requires ln_use_jattr = .false. ' )
+
+         IF(lwp) THEN                     ! Control print
+            WRITE(numout,*)
+            WRITE(numout,*) 'AGRIF boundaries and ghost cells:'
+            WRITE(numout,*) 'lk_west' , lk_west
+            WRITE(numout,*) 'lk_east' , lk_east
+            WRITE(numout,*) 'lk_south', lk_south
+            WRITE(numout,*) 'lk_north', lk_north
+            WRITE(numout,*) 'nbghostcells_y_s', nbghostcells_y_s
+            WRITE(numout,*) 'nbghostcells_y_n', nbghostcells_y_n
+            WRITE(numout,*) 'nbghostcells_x_w', nbghostcells_x_w
+            WRITE(numout,*) 'nbghostcells_x_e', nbghostcells_x_e
+         ENDIF
+      ENDIF
+
+9000  FORMAT (a, i4)
       !
       !
@@ -1106 +1190 @@
           
             bounds_chunks(1,1,1,1) = bounds(1,1,2)
-            bounds_chunks(1,1,2,1) = nn_hls+1
+            bounds_chunks(1,1,2,1) = nn_hls+1 
        
             bounds_chunks(2,1,1,2) = nn_hls+1 
@@ -1161 +1245 @@
       IF( isens == 1 ) THEN
          IF( ptx == 2 ) THEN ! T, V points
-            agrif_external_switch_index = jpiglo-i1+2*nn_hls
+            agrif_external_switch_index = jpiglo-i1+2
          ELSE ! U, F points
-            agrif_external_switch_index = jpiglo-i1+2*nn_hls-1      
+            agrif_external_switch_index = jpiglo-i1+1 
          ENDIF
       ELSE IF( isens ==2 ) THEN

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/BDY/bdydyn2d.F90

@@ -92 +92 @@
                llrecv2(2)   = llrecv2(2)   .OR. lrecv_bdyext(ib_bdy,2,2,ir)     ! might search point towards bdy on the east
                llsend3(3:4) = llsend3(3:4) .OR. lsend_bdyint(ib_bdy,3,3:4,ir)   ! north/south, V points
-               llsend3(3)   = llsend3(3)   .OR. lsend_bdyext(ib_bdy,3,3,ir)     ! neighbour might search point towards its north bdy 
+               llsend3(3)   = llsend3(3)   .OR. lsend_bdyext(ib_bdy,3,3,ir)     ! neighbour might search point towards its north bdy
                llrecv3(3:4) = llrecv3(3:4) .OR. lrecv_bdyint(ib_bdy,3,3:4,ir)   ! north/south, V points
                llrecv3(4)   = llrecv3(4)   .OR. lrecv_bdyext(ib_bdy,3,4,ir)     ! might search point towards bdy on the north

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/BDY/bdyini.F90

@@ -549 +549 @@
                      !
                      icount = icount  + 1
-                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1   ! global to local indexes
-                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1   ! global to local indexes
+                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy) - mig(1) + 1   ! global to local indexes
+                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy) - mjg(1) + 1   ! global to local indexes
                      idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
                      idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
@@ -577 +577 @@
                ! check if point has to be sent     to   a neighbour
                ! W neighbour and on the inner left  side
-               IF( ii == 2     .AND. mpiSnei(nn_hls,jpwe) > -1 )   lsend_bdy(ib_bdy,igrd,jpwe,ir) = .TRUE.
+               IF( ii >= Nis0 .AND. ii < Nis0 + nn_hls .AND. mpiSnei(nn_hls,jpwe) > -1 )   lsend_bdy(ib_bdy,igrd,jpwe,ir) = .TRUE.
                ! E neighbour and on the inner right side
-               IF( ii == jpi-1 .AND. mpiSnei(nn_hls,jpea) > -1 )   lsend_bdy(ib_bdy,igrd,jpea,ir) = .TRUE.
+               IF( ii <= Nie0 .AND. ii > Nie0 - nn_hls .AND. mpiSnei(nn_hls,jpea) > -1 )   lsend_bdy(ib_bdy,igrd,jpea,ir) = .TRUE.
                ! S neighbour and on the inner down side
-               IF( ij == 2     .AND. mpiSnei(nn_hls,jpso) > -1 )   lsend_bdy(ib_bdy,igrd,jpso,ir) = .TRUE.
+               IF( ij >= Njs0 .AND. ij < Njs0 + nn_hls .AND. mpiSnei(nn_hls,jpso) > -1 )   lsend_bdy(ib_bdy,igrd,jpso,ir) = .TRUE.
                ! N neighbour and on the inner up   side
-               IF( ij == jpj-1 .AND. mpiSnei(nn_hls,jpno) > -1 )   lsend_bdy(ib_bdy,igrd,jpno,ir) = .TRUE.
+               IF( ij <= Nje0 .AND. ij > Nje0 - nn_hls .AND. mpiSnei(nn_hls,jpno) > -1 )   lsend_bdy(ib_bdy,igrd,jpno,ir) = .TRUE.
                !
                ! check if point has to be received from a neighbour
                ! W neighbour and on the outter left  side
-               IF( ii == 1     .AND. mpiRnei(nn_hls,jpwe) > -1 )   lrecv_bdy(ib_bdy,igrd,jpwe,ir) = .TRUE.
+               IF( ii  < Nis0                          .AND. mpiRnei(nn_hls,jpwe) > -1 )   lrecv_bdy(ib_bdy,igrd,jpwe,ir) = .TRUE.
                ! E neighbour and on the outter right side
-               IF( ii == jpi   .AND. mpiRnei(nn_hls,jpea) > -1 )   lrecv_bdy(ib_bdy,igrd,jpea,ir) = .TRUE.
+               IF( ii  > Nie0                          .AND. mpiRnei(nn_hls,jpea) > -1 )   lrecv_bdy(ib_bdy,igrd,jpea,ir) = .TRUE.
                ! S neighbour and on the outter down side
-               IF( ij == 1     .AND. mpiRnei(nn_hls,jpso) > -1 )   lrecv_bdy(ib_bdy,igrd,jpso,ir) = .TRUE.
+               IF( ij  < Njs0                          .AND. mpiRnei(nn_hls,jpso) > -1 )   lrecv_bdy(ib_bdy,igrd,jpso,ir) = .TRUE.
                ! N neighbour and on the outter up   side
-               IF( ij == jpj   .AND. mpiRnei(nn_hls,jpno) > -1 )   lrecv_bdy(ib_bdy,igrd,jpno,ir) = .TRUE.
+               IF( ij  > Nje0                          .AND. mpiRnei(nn_hls,jpno) > -1 )   lrecv_bdy(ib_bdy,igrd,jpno,ir) = .TRUE.
                !
             END DO
@@ -642 +642 @@
       ! Read global 2D mask at T-points: bdytmask
       ! -----------------------------------------
-      ! bdytmask = 1  on the computational domain AND on open boundaries
+      ! bdytmask = 1  on the computational domain but not on open boundaries
       !          = 0  elsewhere   
 
@@ -732 +732 @@
                !
                ! search neighbour in the  west/east  direction
-               ! Rim is on the halo and computed ocean is towards exterior of mpi domain  
+               !
+               ! Rim is on the halo and computed ocean is towards exterior of mpi domain :
                !      <--    (o exterior)     -->  
                ! (1)  o|x         OR    (2)   x|o
                !       |___                 ___|
+               ! ==> cannot compute the point x -> need to receive it
                IF( iibi==0     .OR. ii1==0     .OR. ii2==0     .OR. ii3==0     )   lrecv_bdyint(ib_bdy,igrd,jpwe,ir) = .TRUE.
+               IF( iibe==0                                                     )   lrecv_bdyext(ib_bdy,igrd,jpwe,ir) = .TRUE.
                IF( iibi==jpi+1 .OR. ii1==jpi+1 .OR. ii2==jpi+1 .OR. ii3==jpi+1 )   lrecv_bdyint(ib_bdy,igrd,jpea,ir) = .TRUE.  
-               IF( iibe==0                                                     )   lrecv_bdyext(ib_bdy,igrd,jpwe,ir) = .TRUE.
                IF( iibe==jpi+1                                                 )   lrecv_bdyext(ib_bdy,igrd,jpea,ir) = .TRUE.  
-               ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo
+               ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo.
                ! :¨¨¨¨¨|¨¨-->    |                                             |    <--¨¨|¨¨¨¨¨: 
                ! :     |  x:o    |    neighbour limited by ... would need o    |    o:x  |     :
                ! :.....|_._:_____|   (1) W neighbour         E neighbour (2)   |_____:_._|.....:
-               IF( ii==2     .AND. mpiSnei(nn_hls,jpwe) > -1 .AND. &
-                  & ( iibi==3     .OR. ii1==3     .OR. ii2==3     .OR. ii3==3    ) )   lsend_bdyint(ib_bdy,igrd,jpwe,ir) = .TRUE.
-               IF( ii==jpi-1 .AND. mpiSnei(nn_hls,jpea) > -1 .AND. &
-                  & ( iibi==jpi-2 .OR. ii1==jpi-2 .OR. ii2==jpi-2 .OR. ii3==jpi-2) )   lsend_bdyint(ib_bdy,igrd,jpea,ir) = .TRUE.
-               IF( ii==2     .AND. mpiSnei(nn_hls,jpwe) > -1 .AND. iibe==3     )   lsend_bdyext(ib_bdy,igrd,jpwe,ir) = .TRUE.
-               IF( ii==jpi-1 .AND. mpiSnei(nn_hls,jpea) > -1 .AND. iibe==jpi-2 )   lsend_bdyext(ib_bdy,igrd,jpea,ir) = .TRUE.
+               ! ==> the neighbour cannot compute the point x -> need to send it
+               IF( ii ==     2*nn_hls   .AND. mpiSnei(nn_hls,jpwe) > -1 ) THEN   ! 2*nn_hls      -> ji=jpi of western neighbour
+                  IF( iibi==ii+1 .OR. ii1==ii+1 .OR. ii2==ii+1 .OR. ii3==ii+1 )   lsend_bdyint(ib_bdy,igrd,jpwe,ir) = .TRUE.
+                  IF( iibe==ii+1                                              )   lsend_bdyext(ib_bdy,igrd,jpwe,ir) = .TRUE.
+               ENDIF
+               IF( ii == jpi-2*nn_hls+1 .AND. mpiSnei(nn_hls,jpea) > -1 ) THEN   ! jpi-2*nn_hls+1-> ji=1   of eastern neighbour
+                  IF( iibi==ii-1 .OR. ii1==ii-1 .OR. ii2==ii-1 .OR. ii3==ii-1 )   lsend_bdyint(ib_bdy,igrd,jpea,ir) = .TRUE.
+                  IF( iibe==ii-1                                              )   lsend_bdyext(ib_bdy,igrd,jpea,ir) = .TRUE.
+               ENDIF
                !
                ! search neighbour in the north/south direction   
+               !
                ! Rim is on the halo and computed ocean is towards exterior of mpi domain
+               ! ==> cannot compute the point x -> need to receive it
                !(3)   |       |         ^   ___o___     
                !  |   |___x___|   OR    |  |   x   |
                !  v       o           (4)  |       |
                IF( ijbi==0     .OR. ij1==0     .OR. ij2==0     .OR. ij3==0     )   lrecv_bdyint(ib_bdy,igrd,jpso,ir) = .TRUE.
+               IF( ijbe==0                                                     )   lrecv_bdyext(ib_bdy,igrd,jpso,ir) = .TRUE.
                IF( ijbi==jpj+1 .OR. ij1==jpj+1 .OR. ij2==jpj+1 .OR. ij3==jpj+1 )   lrecv_bdyint(ib_bdy,igrd,jpno,ir) = .TRUE.
-               IF( ijbe==0                                                     )   lrecv_bdyext(ib_bdy,igrd,jpso,ir) = .TRUE.
                IF( ijbe==jpj+1                                                 )   lrecv_bdyext(ib_bdy,igrd,jpno,ir) = .TRUE.
                ! Check if neighbour has its rim parallel to its mpi subdomain     _________  border and next to its halo
@@ -764 +771 @@
                !   |  |¨¨¨¨x¨¨¨¨|   neighbour limited by ... would need o     |  |....x....|
                !      :_________:  (3) S neighbour          N neighbour (4)   v  |    o    |   
-               IF( ij==2     .AND. mpiSnei(nn_hls,jpso) > -1 .AND. &
-                  & ( ijbi==3     .OR. ij1==3     .OR. ij2==3     .OR. ij3==3    ) )   lsend_bdyint(ib_bdy,igrd,jpso,ir) = .TRUE.
-               IF( ij==jpj-1 .AND. mpiSnei(nn_hls,jpno) > -1 .AND. &
-                  & ( ijbi==jpj-2 .OR. ij1==jpj-2 .OR. ij2==jpj-2 .OR. ij3==jpj-2) )   lsend_bdyint(ib_bdy,igrd,jpno,ir) = .TRUE.
-               IF( ij==2     .AND. mpiSnei(nn_hls,jpso) > -1 .AND. ijbe==3     )   lsend_bdyext(ib_bdy,igrd,jpso,ir) = .TRUE.
-               IF( ij==jpj-1 .AND. mpiSnei(nn_hls,jpno) > -1 .AND. ijbe==jpj-2 )   lsend_bdyext(ib_bdy,igrd,jpno,ir) = .TRUE.
+               ! ==> the neighbour cannot compute the point x -> need to send it
+               IF( ij ==     2*nn_hls   .AND. mpiSnei(nn_hls,jpso) > -1 ) THEN   ! 2*nn_hls      -> jj=jpj of southern neighbour 
+                  IF( ijbi==ij+1 .OR. ij1==ij+1 .OR. ij2==ij+1 .OR. ij3==ij+1 )   lsend_bdyint(ib_bdy,igrd,jpso,ir) = .TRUE.
+                  IF( ijbe==ij+1                                              )   lsend_bdyext(ib_bdy,igrd,jpso,ir) = .TRUE.
+               ENDIF
+               IF( ij == jpj-2*nn_hls+1 .AND. mpiSnei(nn_hls,jpno) > -1 ) THEN   ! jpj-2*nn_hls+1-> jj=1   of northern neighbour
+                  IF( ijbi==ij-1 .OR. ij1==ij-1 .OR. ij2==ij-1 .OR. ij3==ij-1 )   lsend_bdyint(ib_bdy,igrd,jpno,ir) = .TRUE.
+                  IF( ijbe==ij-1                                              )   lsend_bdyext(ib_bdy,igrd,jpno,ir) = .TRUE.
+               ENDIF
             END DO
          END DO
@@ -813 +823 @@
       !!                - and look at the ocean neighbours to compute ntreat
       !!----------------------------------------------------------------------
-      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pfmask   ! temporary fmask excluding coastal boundary condition (shlat)
-      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pumask, pvmask   ! temporary t/u/v mask array
-      LOGICAL                             , INTENT (in   ) :: lrim0    ! .true. -> rim 0   .false. -> rim 1
+      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pumask, pvmask   ! temporary u/v mask array
+      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pfmask           ! temporary fmask excluding coastal boundary condition (shlat)
+      LOGICAL                             , INTENT (in   ) :: lrim0            ! .true. -> rim 0   .false. -> rim 1
       INTEGER  ::   ib_bdy, ii, ij, igrd, ib, icount       ! dummy loop indices
       INTEGER  ::   i_offset, j_offset, inn                ! local integer
@@ -830 +840 @@
       DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
 
-         ! Calculate relationship of U direction to the local orientation of the boundary
-         ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
-         ! flagu =  0 : u is tangential
-         ! flagu =  1 : u is normal to the boundary and is direction is inward
-         DO igrd = 1, jpbgrd 
-            SELECT CASE( igrd )
-               CASE( 1 )   ;   zmask => pumask     ;   i_offset = 0
-               CASE( 2 )   ;   zmask => bdytmask   ;   i_offset = 1
-               CASE( 3 )   ;   zmask => pfmask     ;   i_offset = 0
-            END SELECT 
-            icount = 0
-            ztmp(:,:) = -999._wp
+         DO igrd = 1, jpbgrd
+            
             IF( lrim0 ) THEN   ! extent of rim 0
                ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
@@ -847 +847 @@
                ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
             END IF
+
+            ! Calculate relationship of U direction to the local orientation of the boundary
+            ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
+            ! flagu =  0 : u is tangential
+            ! flagu =  1 : u is normal to the boundary and is direction is inward
+            SELECT CASE( igrd )
+               CASE( 1 )   ;   zmask => pumask     ;   i_offset = 0   ! U(i-1)   T(i)   U(i  )
+               CASE( 2 )   ;   zmask => bdytmask   ;   i_offset = 1   ! T(i  )   U(i)   T(i+1)
+               CASE( 3 )   ;   zmask => pfmask     ;   i_offset = 0   ! F(i-1)   V(i)   F(i  )
+            END SELECT 
+            icount = 0
+            ztmp(:,:) = -999._wp
             DO ib = ibeg, iend 
                ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE
+               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
                zwfl = zmask(ii+i_offset-1,ij)
                zefl = zmask(ii+i_offset  ,ij)
-               ! This error check only works if you are using the bdyXmask arrays
-               IF( i_offset == 1 .and. zefl + zwfl == 2. ) THEN
+               ! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
+               IF( i_offset == 1 .and. zefl + zwfl == 2._wp ) THEN
                   icount = icount + 1
                   IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
@@ -876 +888 @@
                idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij)
             END DO
-         END DO
-
-         ! Calculate relationship of V direction to the local orientation of the boundary
-         ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
-         ! flagv =  0 : v is tangential
-         ! flagv =  1 : v is normal to the boundary and is direction is inward
-         DO igrd = 1, jpbgrd 
+            
+            ! Calculate relationship of V direction to the local orientation of the boundary
+            ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
+            ! flagv =  0 : v is tangential
+            ! flagv =  1 : v is normal to the boundary and is direction is inward
             SELECT CASE( igrd )
                CASE( 1 )   ;   zmask => pvmask     ;   j_offset = 0
@@ -890 +900 @@
             icount = 0
             ztmp(:,:) = -999._wp
-            IF( lrim0 ) THEN   ! extent of rim 0
-               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
-            ELSE               ! extent of rim 1
-               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
-            END IF
             DO ib = ibeg, iend
                ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE
+               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
                zsfl = zmask(ii,ij+j_offset-1)
                znfl = zmask(ii,ij+j_offset  )
-               ! This error check only works if you are using the bdyXmask arrays
-               IF( j_offset == 1 .and. znfl + zsfl == 2. ) THEN
+               ! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
+               IF( j_offset == 1 .and. znfl + zsfl == 2._wp ) THEN
                   IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
                   icount = icount + 1
@@ -924 +929 @@
                idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij)
             END DO
-         END DO
-         !
-      END DO ! ib_bdy
       
-      DO ib_bdy = 1, nb_bdy
-         DO igrd = 1, jpbgrd
+            ! Calculate ntreat
             SELECT CASE( igrd )
                CASE( 1 )   ;   zmask => bdytmask 
@@ -936 +937 @@
             END SELECT
             ztmp(:,:) = -999._wp
-            IF( lrim0 ) THEN   ! extent of rim 0
-               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
-            ELSE               ! extent of rim 1
-               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
-            END IF
             DO ib = ibeg, iend
                ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE
-               llnon = zmask(ii  ,ij+1) == 1.  
-               llson = zmask(ii  ,ij-1) == 1. 
-               llean = zmask(ii+1,ij  ) == 1. 
-               llwen = zmask(ii-1,ij  ) == 1. 
+               IF( ii < Nis0 .OR. ii > Nie0 .OR. ij < Njs0 .OR. ij > Nje0 )  CYCLE   ! call lbc_lnk -> no need to compute these pts
+               llnon = zmask(ii  ,ij+1) == 1._wp  
+               llson = zmask(ii  ,ij-1) == 1._wp 
+               llean = zmask(ii+1,ij  ) == 1._wp 
+               llwen = zmask(ii-1,ij  ) == 1._wp 
                inn  = COUNT( (/ llnon, llson, llean, llwen /) )
                IF( inn == 0 ) THEN   ! no neighbours -> interior of a corner  or  cluster of rim points
@@ -954 +950 @@
                   !  1 |   o      !  2  o   |    !  3 | x        !  4     x |   !      |   |   -> error
                   !    |_x_ _     !    _ _x_|    !    |   o      !      o   |   !      |x_x|
-                  IF(     zmask(ii+1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 1.
-                  ELSEIF( zmask(ii-1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 2.
-                  ELSEIF( zmask(ii+1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 3.
-                  ELSEIF( zmask(ii-1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 4.
-                  ELSE                                    ;   ztmp(ii,ij) = -1.
+                  IF(     zmask(ii+1,ij+1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 1._wp
+                  ELSEIF( zmask(ii-1,ij+1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 2._wp
+                  ELSEIF( zmask(ii+1,ij-1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 3._wp
+                  ELSEIF( zmask(ii-1,ij-1) == 1._wp ) THEN   ;   ztmp(ii,ij) = 4._wp
+                  ELSE                                       ;   ztmp(ii,ij) = -1._wp
                      WRITE(ctmp1,*) 'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
                        ' on boundary set ', ib_bdy, ' has no free ocean neighbour'
@@ -973 +969 @@
                   ! 5  | x o     ! 6   o x |   ! 7  __x__   ! 8    x
                   !    |         !         |   !            !      o
-                  IF( llean )   ztmp(ii,ij) = 5.
-                  IF( llwen )   ztmp(ii,ij) = 6.
-                  IF( llnon )   ztmp(ii,ij) = 7.
-                  IF( llson )   ztmp(ii,ij) = 8.
+                  IF( llean )   ztmp(ii,ij) = 5._wp
+                  IF( llwen )   ztmp(ii,ij) = 6._wp
+                  IF( llnon )   ztmp(ii,ij) = 7._wp
+                  IF( llson )   ztmp(ii,ij) = 8._wp
                END IF
                IF( inn == 2 ) THEN   ! exterior of a corner
@@ -982 +978 @@
                   !  9 ____x o    ! 10   o x___   ! 11     x o      ! 12   o x      
                   !         |     !       |       !        o        !        o 
-                  IF( llnon .AND. llean )   ztmp(ii,ij) =  9.
-                  IF( llnon .AND. llwen )   ztmp(ii,ij) = 10.
-                  IF( llson .AND. llean )   ztmp(ii,ij) = 11.
-                  IF( llson .AND. llwen )   ztmp(ii,ij) = 12.
+                  IF( llnon .AND. llean )   ztmp(ii,ij) =  9._wp
+                  IF( llnon .AND. llwen )   ztmp(ii,ij) = 10._wp
+                  IF( llson .AND. llean )   ztmp(ii,ij) = 11._wp
+                  IF( llson .AND. llwen )   ztmp(ii,ij) = 12._wp
                END IF
                IF( inn == 3 ) THEN   ! 3 neighbours     __   __
@@ -991 +987 @@
                   ! 13  _| x o    ! 14   o x |_   ! 15   o x o    ! 16  o x o       
                   !    |   o      !        o   |  !        o      !    __|¨|__    
-                  IF( llnon .AND. llean .AND. llson )   ztmp(ii,ij) = 13.
-                  IF( llnon .AND. llwen .AND. llson )   ztmp(ii,ij) = 14.
-                  IF( llwen .AND. llson .AND. llean )   ztmp(ii,ij) = 15.
-                  IF( llwen .AND. llnon .AND. llean )   ztmp(ii,ij) = 16.
+                  IF( llnon .AND. llean .AND. llson )   ztmp(ii,ij) = 13._wp
+                  IF( llnon .AND. llwen .AND. llson )   ztmp(ii,ij) = 14._wp
+                  IF( llwen .AND. llson .AND. llean )   ztmp(ii,ij) = 15._wp
+                  IF( llwen .AND. llnon .AND. llean )   ztmp(ii,ij) = 16._wp
                END IF
                IF( inn == 4 ) THEN
@@ -1012 +1008 @@
                idx_bdy(ib_bdy)%ntreat(ib,igrd) = NINT(ztmp(ii,ij))
             END DO
-         END DO
-      END DO
+            !
+         END DO   ! jpbgrd
+         !
+      END DO   ! ib_bdy
 
     END SUBROUTINE bdy_rim_treat

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/DIA/dia25h.F90

@@ -32 +32 @@
    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   en_25h  , rmxln_25h
 
+!! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -51 +53 @@
       INTEGER ::   ios                 ! Local integer output status for namelist read
       INTEGER ::   ierror              ! Local integer for memory allocation
+      INTEGER ::   ji, jj, jk
       !
       NAMELIST/nam_dia25h/ ln_dia25h
@@ -73 +76 @@
       ! ------------------- !
       !                                ! ocean arrays
-      ALLOCATE( tn_25h (jpi,jpj,jpk), sn_25h (jpi,jpj,jpk), sshn_25h(jpi,jpj)  ,     &
-         &      un_25h (jpi,jpj,jpk), vn_25h (jpi,jpj,jpk), wn_25h(jpi,jpj,jpk),     &
-         &      avt_25h(jpi,jpj,jpk), avm_25h(jpi,jpj,jpk),                      STAT=ierror )
+      ALLOCATE( tn_25h (A2D(0),jpk), sn_25h (A2D(0),jpk), sshn_25h(A2D(0))  ,     &
+         &      un_25h (A2D(0),jpk), vn_25h (A2D(0),jpk), wn_25h(A2D(0),jpk),     &
+         &      avt_25h(A2D(0),jpk), avm_25h(A2D(0),jpk),                      STAT=ierror )
       IF( ierror > 0 ) THEN
          CALL ctl_stop( 'dia_25h: unable to allocate ocean arrays' )   ;   RETURN
       ENDIF
       IF( ln_zdftke ) THEN             ! TKE physics
-         ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror )
+         ALLOCATE( en_25h(A2D(0),jpk), STAT=ierror )
          IF( ierror > 0 ) THEN
             CALL ctl_stop( 'dia_25h: unable to allocate en_25h' )   ;   RETURN
@@ -86 +89 @@
       ENDIF
       IF( ln_zdfgls ) THEN             ! GLS physics
-         ALLOCATE( en_25h(jpi,jpj,jpk), rmxln_25h(jpi,jpj,jpk), STAT=ierror )
+         ALLOCATE( en_25h(A2D(0),jpk), rmxln_25h(A2D(0),jpk), STAT=ierror )
          IF( ierror > 0 ) THEN
             CALL ctl_stop( 'dia_25h: unable to allocate en_25h and rmxln_25h' )   ;   RETURN
@@ -94 +97 @@
       ! 2 - Assign Initial Values !
       ! ------------------------- !
-      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 
-      tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kbb)
-      sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kbb)
-      sshn_25h(:,:)   = ssh(:,:,Kbb)
-      un_25h  (:,:,:) = uu  (:,:,:,Kbb)
-      vn_25h  (:,:,:) = vv  (:,:,:,Kbb)
-      avt_25h (:,:,:) = avt (:,:,:)
-      avm_25h (:,:,:) = avm (:,:,:)
+      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible)
+      DO_3D( 0, 0, 0, 0, 1, jpk )
+         tn_25h (ji,jj,jk) = ts (ji,jj,jk,jp_tem,Kbb)
+         sn_25h (ji,jj,jk) = ts (ji,jj,jk,jp_sal,Kbb)
+         un_25h (ji,jj,jk) = uu (ji,jj,jk,Kbb)
+         vn_25h (ji,jj,jk) = vv (ji,jj,jk,Kbb)
+         avt_25h(ji,jj,jk) = avt(ji,jj,jk)
+         avm_25h(ji,jj,jk) = avm(ji,jj,jk)
+      END_3D
+      DO_2D( 0, 0, 0, 0 )
+         sshn_25h(ji,jj) = ssh(ji,jj,Kbb)
+      END_2D
       IF( ln_zdftke ) THEN
-         en_25h(:,:,:) = en(:,:,:)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            en_25h(ji,jj,jk) = en(ji,jj,jk)
+         END_3D
       ENDIF
       IF( ln_zdfgls ) THEN
-         en_25h   (:,:,:) = en    (:,:,:)
-         rmxln_25h(:,:,:) = hmxl_n(:,:,:)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            en_25h   (ji,jj,jk) = en    (ji,jj,jk)
+            rmxln_25h(ji,jj,jk) = hmxl_n(ji,jj,jk)
+         END_3D
       ENDIF
 #if defined key_si3
@@ -132 +143 @@
       REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi   ! local scalars
       INTEGER                          ::   i_steps                           ! no of timesteps per hour
-      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                ! workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                              ! workspace
-      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                             ! workspace
+      REAL(wp), DIMENSION(A2D(0)    )  ::   zw2d, un_dm, vn_dm                ! workspace
+      REAL(wp), DIMENSION(A2D(0),jpk)  ::   zw3d                              ! workspace
+      REAL(wp), DIMENSION(A2D(0),3)    ::   zwtmb                             ! workspace
       !!----------------------------------------------------------------------
 
@@ -151 +162 @@
       ! wn_25h could not be initialised in dia_25h_init, so we do it here instead
       IF( kt == nn_it000 ) THEN
-         wn_25h(:,:,:) = ww(:,:,:)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            wn_25h(ji,jj,jk) = ww(ji,jj,jk)
+         END_3D
       ENDIF
 
@@ -162 +175 @@
          ENDIF
 
-         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + ts (:,:,:,jp_tem,Kmm)
-         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + ts (:,:,:,jp_sal,Kmm)
-         sshn_25h(:,:)       = sshn_25h(:,:)   + ssh(:,:,Kmm)
-         un_25h  (:,:,:)     = un_25h  (:,:,:) + uu  (:,:,:,Kmm)
-         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vv  (:,:,:,Kmm)
-         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + ww  (:,:,:)
-         avt_25h (:,:,:)     = avt_25h (:,:,:) + avt (:,:,:)
-         avm_25h (:,:,:)     = avm_25h (:,:,:) + avm (:,:,:)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            tn_25h  (ji,jj,jk) = tn_25h  (ji,jj,jk) + ts (ji,jj,jk,jp_tem,Kmm)
+            sn_25h  (ji,jj,jk) = sn_25h  (ji,jj,jk) + ts (ji,jj,jk,jp_sal,Kmm)
+            un_25h  (ji,jj,jk) = un_25h  (ji,jj,jk) + uu (ji,jj,jk,Kmm)
+            vn_25h  (ji,jj,jk) = vn_25h  (ji,jj,jk) + vv (ji,jj,jk,Kmm)
+            wn_25h  (ji,jj,jk) = wn_25h  (ji,jj,jk) + ww (ji,jj,jk)
+            avt_25h (ji,jj,jk) = avt_25h (ji,jj,jk) + avt(ji,jj,jk)
+            avm_25h (ji,jj,jk) = avm_25h (ji,jj,jk) + avm(ji,jj,jk)
+         END_3D
+         DO_2D( 0, 0, 0, 0 )
+            sshn_25h(ji,jj)    = sshn_25h(ji,jj)    + ssh(ji,jj,Kmm)
+         END_2D
          IF( ln_zdftke ) THEN
-            en_25h(:,:,:)    = en_25h  (:,:,:) + en(:,:,:)
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               en_25h(ji,jj,jk) = en_25h(ji,jj,jk) + en(ji,jj,jk)
+            END_3D
          ENDIF
          IF( ln_zdfgls ) THEN
-            en_25h   (:,:,:) = en_25h   (:,:,:) + en    (:,:,:)
-            rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + hmxl_n(:,:,:)
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               en_25h   (ji,jj,jk) = en_25h   (ji,jj,jk) + en    (ji,jj,jk)
+               rmxln_25h(ji,jj,jk) = rmxln_25h(ji,jj,jk) + hmxl_n(ji,jj,jk)
+            END_3D
          ENDIF
          cnt_25h = cnt_25h + 1
@@ -212 +233 @@
          zmdi=1.e+20 !missing data indicator for masking
          ! write tracers (instantaneous)
-         zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = tn_25h(ji,jj,jk)*tmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+         END_3D
          CALL iom_put("temper25h", zw3d)   ! potential temperature
-         zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = sn_25h(ji,jj,jk)*tmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+         END_3D
          CALL iom_put( "salin25h", zw3d  )   ! salinity
-         zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+         DO_2D( 0, 0, 0, 0 )
+            zw2d(ji,jj) = sshn_25h(ji,jj)*tmask(ji,jj,1) + zmdi*(1.0-tmask(ji,jj,1))
+         END_2D
          IF( ll_wd ) THEN
             CALL iom_put( "ssh25h", zw2d+ssh_ref )   ! sea surface 
@@ -223 +250 @@
          ENDIF
          ! Write velocities (instantaneous)
-         zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = un_25h(ji,jj,jk)*umask(ji,jj,jk) + zmdi*(1.0-umask(ji,jj,jk))
+         END_3D
          CALL iom_put("vozocrtx25h", zw3d)    ! i-current
-         zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = vn_25h(ji,jj,jk)*vmask(ji,jj,jk) + zmdi*(1.0-vmask(ji,jj,jk))
+         END_3D
          CALL iom_put("vomecrty25h", zw3d  )   ! j-current
-         zw3d(:,:,:) = wn_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = wn_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+         END_3D
          CALL iom_put("vovecrtz25h", zw3d )   ! k-current
          ! Write vertical physics
-         zw3d(:,:,:) = avt_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = avt_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+         END_3D
          CALL iom_put("avt25h", zw3d )   ! diffusivity
-         zw3d(:,:,:) = avm_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = avm_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+         END_3D
          CALL iom_put("avm25h", zw3d)   ! viscosity
          IF( ln_zdftke ) THEN
-            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               zw3d(ji,jj,jk) = en_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+            END_3D
             CALL iom_put("tke25h", zw3d)   ! tke
          ENDIF
          IF( ln_zdfgls ) THEN
-            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               zw3d(ji,jj,jk) = en_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+            END_3D
             CALL iom_put("tke25h", zw3d)   ! tke
-            zw3d(:,:,:) = rmxln_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               zw3d(ji,jj,jk) = rmxln_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+            END_3D
             CALL iom_put( "mxln25h",zw3d)
          ENDIF
          !
          ! After the write reset the values to cnt=1 and sum values equal current value 
-         tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kmm)
-         sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kmm)
-         sshn_25h(:,:)   = ssh(:,:,Kmm)
-         un_25h  (:,:,:) = uu  (:,:,:,Kmm)
-         vn_25h  (:,:,:) = vv  (:,:,:,Kmm)
-         wn_25h  (:,:,:) = ww  (:,:,:)
-         avt_25h (:,:,:) = avt (:,:,:)
-         avm_25h (:,:,:) = avm (:,:,:)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            tn_25h  (ji,jj,jk) = ts (ji,jj,jk,jp_tem,Kmm)
+            sn_25h  (ji,jj,jk) = ts (ji,jj,jk,jp_sal,Kmm)
+            un_25h  (ji,jj,jk) = uu (ji,jj,jk,Kmm)
+            vn_25h  (ji,jj,jk) = vv (ji,jj,jk,Kmm)
+            wn_25h  (ji,jj,jk) = ww (ji,jj,jk)
+            avt_25h (ji,jj,jk) = avt(ji,jj,jk)
+            avm_25h (ji,jj,jk) = avm(ji,jj,jk)
+         END_3D
+         DO_2D( 0, 0, 0, 0 )
+            sshn_25h(ji,jj)    = ssh(ji,jj,Kmm)
+         END_2D
          IF( ln_zdftke ) THEN
-            en_25h(:,:,:) = en(:,:,:)
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               en_25h(ji,jj,jk) = en(ji,jj,jk)
+            END_3D
          ENDIF
          IF( ln_zdfgls ) THEN
-            en_25h   (:,:,:) = en    (:,:,:)
-            rmxln_25h(:,:,:) = hmxl_n(:,:,:)
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               en_25h   (ji,jj,jk) = en    (ji,jj,jk)
+               rmxln_25h(ji,jj,jk) = hmxl_n(ji,jj,jk)
+            END_3D
          ENDIF
          cnt_25h = 1

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/DIA/diahth.F90

@@ -334 +334 @@
    SUBROUTINE dia_hth_htc( Kmm, pdep, pt, phtc )
       !
-      INTEGER , INTENT(in) :: Kmm      ! ocean time level index
-      REAL(wp), INTENT(in) :: pdep     ! depth over the heat content
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pt   
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc  
-      !
-      INTEGER  :: ji, jj, jk, ik
-      REAL(wp), DIMENSION(jpi,jpj) :: zthick
-      INTEGER , DIMENSION(jpi,jpj) :: ilevel
+      INTEGER , INTENT(in) ::   Kmm      ! ocean time level index
+      REAL(wp), INTENT(in) ::   pdep     ! depth over the heat content
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)    ::   pt   
+      REAL(wp), DIMENSION(jpi,jpj),     INTENT(inout) ::   phtc  
+      !
+      INTEGER  ::   ji, jj, jk, ik
+      REAL(wp), DIMENSION(jpi,jpj) ::   zthick
+      INTEGER , DIMENSION(jpi,jpj) ::   ilevel
 
 
       ! surface boundary condition
       
-      IF( .NOT. ln_linssh ) THEN   ;   zthick(:,:) = 0._wp       ;   phtc(:,:) = 0._wp                                   
+      IF( .NOT. ln_linssh ) THEN   ;   zthick(:,:) = 0._wp          ;   phtc(:,:) = 0._wp                                   
       ELSE                         ;   zthick(:,:) = ssh(:,:,Kmm)   ;   phtc(:,:) = pt(:,:,1) * ssh(:,:,Kmm) * tmask(:,:,1)   
       ENDIF
       !
       ilevel(:,:) = 1
-      DO_3D( 1, 1, 1, 1, 2, jpkm1 )
-         IF( ( gdept(ji,jj,jk,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
-             ilevel(ji,jj) = jk
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+         IF( ( gdepw(ji,jj,jk+1,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
+             ilevel(ji,jj) = jk+1
              zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm)
              phtc  (ji,jj) = phtc  (ji,jj) + e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk)
@@ -361 +361 @@
       DO_2D( 1, 1, 1, 1 )
          ik = ilevel(ji,jj)
-         zthick(ji,jj) = pdep - zthick(ji,jj)   !   remaining thickness to reach depht pdep
-         phtc(ji,jj)   = phtc(ji,jj)    &
-            &           + pt (ji,jj,ik+1) * MIN( e3t(ji,jj,ik+1,Kmm), zthick(ji,jj) ) &
-                                                       * tmask(ji,jj,ik+1)
+         IF( tmask(ji,jj,ik) == 1 ) THEN
+            zthick(ji,jj) = MIN ( gdepw(ji,jj,ik+1,Kmm), pdep ) - zthick(ji,jj)   ! remaining thickness to reach dephw pdep
+            phtc(ji,jj)   = phtc(ji,jj) + pt(ji,jj,ik) * zthick(ji,jj)
+         ENDIF
       END_2D
-      !
       !
    END SUBROUTINE dia_hth_htc

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/DIA/diawri.F90

@@ -583 +583 @@
       REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
       !
-      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d, ze3t, zgdept       ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj    ) :: z2d     ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d     ! 3D workspace
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace
       !!----------------------------------------------------------------------
@@ -624 +624 @@
       it = kt
       itmod = kt - nit000 + 1
-
-      ! store e3t for subsitute
-      DO jk = 1, jpk
-         ze3t  (:,:,jk) =  e3t  (:,:,jk,Kmm)
-         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm)
-      END DO
-
 
       ! 1. Define NETCDF files and fields at beginning of first time step
@@ -786 +779 @@
          CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         IF( ALLOCATED(hmld) ) THEN   ! zdf_mxl not called by SWE
+            CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+            CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
+               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+         ENDIF
          CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
@@ -942 +937 @@
 
       IF( .NOT.ln_linssh ) THEN
-         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! heat content
-         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! salt content
-         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
-         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            z3d(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm)
+         END_3D
+         CALL histwrite( nid_T, "votemper", it, z3d, ndim_T , ndex_T  )   ! heat content
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            z3d(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm)
+         END_3D
+         CALL histwrite( nid_T, "vosaline", it, z3d, ndim_T , ndex_T  )   ! salt content
+         DO_2D( 0, 0, 0, 0 )
+            z2d(ji,jj   ) = ts(ji,jj, 1,jp_tem,Kmm) * e3t(ji,jj, 1,Kmm)
+         END_2D
+         CALL histwrite( nid_T, "sosstsst", it, z2d, ndim_hT, ndex_hT )   ! sea surface heat content
+         DO_2D( 0, 0, 0, 0 )
+            z2d(ji,jj   ) = ts(ji,jj, 1,jp_sal,Kmm) * e3t(ji,jj, 1,Kmm)
+         END_2D
+         CALL histwrite( nid_T, "sosaline", it, z2d, ndim_hT, ndex_hT )   ! sea surface salinity content
       ELSE
          CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature
@@ -953 +960 @@
       ENDIF
       IF( .NOT.ln_linssh ) THEN
-         zw3d(:,:,:) = ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
-         CALL histwrite( nid_T, "vovvle3t", it, ze3t (:,:,:)     , ndim_T , ndex_T  )   ! level thickness
-         CALL histwrite( nid_T, "vovvldep", it, zgdept , ndim_T , ndex_T  )   ! t-point depth 
-         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
-      ENDIF
-      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)          , ndim_hT, ndex_hT )   ! sea surface height
-      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+           z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm)     ! 3D workspace for qco substitution
+         END_3D
+         CALL histwrite( nid_T, "vovvle3t", it, z3d        , ndim_T , ndex_T  )   ! level thickness
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+           z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm)   ! 3D workspace for qco substitution
+         END_3D
+         CALL histwrite( nid_T, "vovvldep", it, z3d        , ndim_T , ndex_T  )   ! t-point depth 
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            z3d(ji,jj,jk) = ( ( e3t(ji,jj,jk,Kmm) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
+         END_3D         
+         CALL histwrite( nid_T, "vovvldef", it, z3d        , ndim_T , ndex_T  )   ! level thickness deformation
+      ENDIF
+      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)  , ndim_hT, ndex_hT )   ! sea surface height
+      DO_2D( 0, 0, 0, 0 )
+         z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
+      END_2D
+      CALL histwrite( nid_T, "sowaflup", it, z2d           , ndim_hT, ndex_hT )   ! upward water flux
       CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
       CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux 
@@ -965 +983 @@
                                                                                   ! in linear free surface case)
       IF( ln_linssh ) THEN
-         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm)
-         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
-         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm)
-         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
-      ENDIF
-      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
+         DO_2D( 0, 0, 0, 0 )
+            z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
+         END_2D
+         CALL histwrite( nid_T, "sosst_cd", it, z2d, ndim_hT, ndex_hT )          ! c/d term on sst
+         DO_2D( 0, 0, 0, 0 )
+            z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
+         END_2D
+         CALL histwrite( nid_T, "sosss_cd", it, z2d, ndim_hT, ndex_hT )          ! c/d term on sss
+      ENDIF
+      DO_2D( 0, 0, 0, 0 )
+         z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
+      END_2D
+      CALL histwrite( nid_T, "sohefldo", it, z2d           , ndim_hT, ndex_hT )   ! total heat flux
       CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
-      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
-      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
+      IF( ALLOCATED(hmld) ) THEN   ! zdf_mxl not called by SWE
+         CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
+         CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
+      ENDIF
       CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
       CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
@@ -1026 +1053 @@
          CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
          CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
-         zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
-         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
+         DO_2D( 0, 0, 0, 0 )
+            z2d(ji,jj) = erp(ji,jj) * ts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1)
+         END_2D
+         CALL histwrite( nid_T, "sosafldp", it, z2d           , ndim_hT, ndex_hT )   ! salt flux damping
       ENDIF
 !      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
@@ -1039 +1068 @@
 #endif
 
-      CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm)            , ndim_U , ndex_U )    ! i-current
+      CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm) , ndim_U , ndex_U )    ! i-current
       CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
 
-      CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm)            , ndim_V , ndex_V  )   ! j-current
+      CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm) , ndim_V , ndex_V  )   ! j-current
       CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
 
       IF( ln_zad_Aimp ) THEN
-         CALL histwrite( nid_W, "vovecrtz", it, ww + wi     , ndim_T, ndex_T )    ! vert. current
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+           z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
+         END_3D         
+         CALL histwrite( nid_W, "vovecrtz", it, z3d         , ndim_T, ndex_T )    ! vert. current
       ELSE
          CALL histwrite( nid_W, "vovecrtz", it, ww          , ndim_T, ndex_T )    ! vert. current
@@ -1093 +1125 @@
       CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
       !!
-      INTEGER :: inum, jk
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept       ! 3D workspace for qco substitution
+      INTEGER ::   ji, jj, jk       ! dummy loop indices
+      INTEGER ::   inum
+      REAL(wp), DIMENSION(jpi,jpj)     :: z2d      
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d      
       !!----------------------------------------------------------------------
       ! 
@@ -1104 +1138 @@
       ENDIF 
       !
-      DO jk = 1, jpk
-         ze3t(:,:,jk) =  e3t(:,:,jk,Kmm)
-         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm)
-      END DO
-      !
       CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
       !
       CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature
       CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity
-      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height
-      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity
-      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity
+      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)         )    ! sea surface height
+      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)        )    ! now i-velocity
+      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)        )    ! now j-velocity
       IF( ln_zad_Aimp ) THEN
-         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi        )    ! now k-velocity
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+           z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
+         END_3D
+         CALL iom_rstput( 0, 0, inum, 'vovecrtz', z3d            )    ! now k-velocity
       ELSE
          CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity
       ENDIF
-      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity
+      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )
       CALL iom_rstput( 0, 0, inum, 'ht'     , ht(:,:)            )    ! now water column height
       !
       IF ( ln_isf ) THEN
          IF (ln_isfcav_mlt) THEN
-            CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )
+            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )
+            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )
+            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) )
+            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) )
             CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 )
          END IF
          IF (ln_isfpar_mlt) THEN
-            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) )    ! now k-velocity
-            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) )
+            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )
+            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )
+            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )
+            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) )
+            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) )
             CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 )
          END IF
@@ -1152 +1184 @@
          CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
       ENDIF
-      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
-      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
+      DO_2D( 0, 0, 0, 0 )
+         z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
+      END_2D
+      CALL iom_rstput( 0, 0, inum, 'sowaflup', z2d               )    ! freshwater budget
+      DO_2D( 0, 0, 0, 0 )
+         z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
+      END_2D
+      CALL iom_rstput( 0, 0, inum, 'sohefldo', z2d               )    ! total heat flux
       CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
       CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
       CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
       CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
-      IF(  .NOT.ln_linssh  ) THEN             
-         CALL iom_rstput( 0, 0, inum, 'vovvldep', zgdept        )    !  T-cell depth 
-         CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t          )    !  T-cell thickness  
+      IF(  .NOT.ln_linssh  ) THEN
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+           z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm)   ! 3D workspace for qco substitution
+         END_3D
+         CALL iom_rstput( 0, 0, inum, 'vovvldep', z3d            )    !  T-cell depth 
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+           z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm)     ! 3D workspace for qco substitution
+         END_3D
+         CALL iom_rstput( 0, 0, inum, 'vovvle3t', z3d            )    !  T-cell thickness  
       END IF
       IF( ln_wave .AND. ln_sdw ) THEN

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/DOM/dom_oce.F90

@@ -274 +274 @@
       !
       ii = ii+1
-      ALLOCATE( mig(jpi), mjg(jpj), mig0(jpi), mjg0(jpj), STAT=ierr(ii) )
-         !
-      ii = ii+1
-      ALLOCATE( mi0(jpiglo), mi1(jpiglo), mj0(jpjglo), mj1(jpjglo), STAT=ierr(ii) )
-         !
-      ii = ii+1
       ALLOCATE( glamt(jpi,jpj) ,    glamu(jpi,jpj) ,  glamv(jpi,jpj) ,  glamf(jpi,jpj) ,     &
          &      gphit(jpi,jpj) ,    gphiu(jpi,jpj) ,  gphiv(jpi,jpj) ,  gphif(jpi,jpj) ,     &

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/DOM/domain.F90

@@ -20 +20 @@
    !!----------------------------------------------------------------------
    !!   dom_init      : initialize the space and time domain
-   !!   dom_glo       : initialize global domain <--> local domain indices
    !!   dom_nam       : read and contral domain namelists
    !!   dom_ctl       : control print for the ocean domain
@@ -124 +123 @@
       !           !==  Reference coordinate system  ==!
       !
-      CALL dom_glo                      ! global domain versus local domain
       CALL dom_nam                      ! read namelist ( namrun, namdom )
       CALL dom_tile_init                ! Tile domain
@@ -244 +242 @@
       !
    END SUBROUTINE dom_init
-
-
-   SUBROUTINE dom_glo
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE dom_glo  ***
-      !!
-      !! ** Purpose :   initialization of global domain <--> local domain indices
-      !!
-      !! ** Method  :
-      !!
-      !! ** Action  : - mig , mjg : local  domain indices ==> global domain, including halos, indices
-      !!              - mig0, mjg0: local  domain indices ==> global domain, excluding halos, indices
-      !!              - mi0 , mi1 : global domain indices ==> local  domain indices
-      !!              - mj0 , mj1   (if global point not in the local domain ==> mi0>mi1 and/or mj0>mj1)
-      !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj   ! dummy loop argument
-      !!----------------------------------------------------------------------
-      !
-      DO ji = 1, jpi                 ! local domain indices ==> global domain indices, including halos
-        mig(ji) = ji + nimpp - 1
-      END DO
-      DO jj = 1, jpj
-        mjg(jj) = jj + njmpp - 1
-      END DO
-      !                              ! local domain indices ==> global domain indices, excluding halos
-      !
-      mig0(:) = mig(:) - nn_hls
-      mjg0(:) = mjg(:) - nn_hls
-      !                              ! global domain, including halos, indices ==> local domain indices
-      !                                   ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the
-      !                                   ! local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain.
-      DO ji = 1, jpiglo
-        mi0(ji) = MAX( 1 , MIN( ji - nimpp + 1, jpi+1 ) )
-        mi1(ji) = MAX( 0 , MIN( ji - nimpp + 1, jpi   ) )
-      END DO
-      DO jj = 1, jpjglo
-        mj0(jj) = MAX( 1 , MIN( jj - njmpp + 1, jpj+1 ) )
-        mj1(jj) = MAX( 0 , MIN( jj - njmpp + 1, jpj   ) )
-      END DO
-      IF(lwp) THEN                   ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_glo : domain: global <<==>> local '
-         WRITE(numout,*) '~~~~~~~ '
-         WRITE(numout,*) '   global domain:   jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpkglo = ', jpkglo
-         WRITE(numout,*) '   local  domain:   jpi    = ', jpi   , ' jpj    = ', jpj   , ' jpk    = ', jpk
-         WRITE(numout,*)
-      ENDIF
-      !
-   END SUBROUTINE dom_glo
 
 

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/DOM/dommsk.F90

@@ -73 +73 @@
       !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
       !!
-      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
-      !!                rows/lines due to cyclic or North Fold boundaries as well
-      !!                as MPP halos.
-      !!      tmask_h : halo mask at t-point, i.e. excluding duplicated rows/lines
-      !!                due to cyclic or North Fold boundaries as well as MPP halos.
-      !!
       !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask 
       !!                         at t-, u-, v- w, wu-, and wv-points (=0. or 1.)
       !!               fmask   : land/ocean mask at f-point (=0., or =1., or 
       !!                         =rn_shlat along lateral boundaries)
-      !!               tmask_i : interior ocean mask 
-      !!               tmask_h : halo mask
-      !!               ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask
+      !!               ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask, i.e. at least 1 wet cell in the vertical
+      !!               tmask_h : halo mask at t-point, i.e. excluding all duplicated rows/lines
+      !!                         due to cyclic or North Fold boundaries as well as MPP halos.
+      !!               tmask_i : ssmask * tmask_h
       !!----------------------------------------------------------------------
       INTEGER, DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! first and last ocean level

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/DYN/dynldf_lap_blp_lf.F90

@@ -102 +102 @@
          DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )                           ! Horizontal slab
             !                                      ! ahm * e3 * curl  (warning: computed for ji-1,jj-1)
-            zcur     = ahmf(ji,jj,jk) * e3f(ji,jj,jk) * r1_e1e2f(ji,jj)               &   ! ahmf already * by fmask   
-               &       * ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
-               &       - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) )
-            zcur_jm1 = ahmf(ji  ,jj-1,jk) * e3f(ji,jj-1,jk) * r1_e1e2f(ji,jj-1)               &   ! ahmf already * by fmask
+            zcur     = ahmf(ji  ,jj  ,jk) * e3f(ji  ,jj  ,jk) * r1_e1e2f(ji  ,jj  )               &   ! ahmf already * by fmask   
+               &       * ( e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
+               &         - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) )
+            zcur_jm1 = ahmf(ji  ,jj-1,jk) * e3f(ji  ,jj-1,jk) * r1_e1e2f(ji  ,jj-1)               &   ! ahmf already * by fmask
                &       * ( e2v(ji+1,jj-1) * pv(ji+1,jj-1,jk) - e2v(ji,jj-1) * pv(ji,jj-1,jk)  &
-               &       - e1u(ji,jj) * pu(ji,jj,jk) + e1u(ji,jj-1) * pu(ji,jj-1,jk) )
-            zcur_im1 = ahmf(ji-1,jj,jk) * e3f(ji-1,jj,jk) * r1_e1e2f(ji-1,jj)         &   ! ahmf already * by fmask
-               &       * ( e2v(ji,jj) * pv(ji,jj,jk) - e2v(ji-1,jj) * pv(ji-1,jj,jk)  &
-               &       - e1u(ji-1,jj+1) * pu(ji-1,jj+1,jk) + e1u(ji-1,jj) * pu(ji-1,jj,jk) )
+               &         - e1u(ji  ,jj  ) * pu(ji  ,jj  ,jk) + e1u(ji,jj-1) * pu(ji,jj-1,jk) )
+            zcur_im1 = ahmf(ji-1,jj  ,jk) * e3f(ji-1,jj  ,jk) * r1_e1e2f(ji-1,jj  )         &   ! ahmf already * by fmask
+               &       * ( e2v(ji  ,jj  ) * pv(ji  ,jj  ,jk) - e2v(ji-1,jj) * pv(ji-1,jj,jk)  &
+               &         - e1u(ji-1,jj+1) * pu(ji-1,jj+1,jk) + e1u(ji-1,jj) * pu(ji-1,jj,jk) )
             !                                      ! ahm * div        (warning: computed for ji,jj)
             zdiv     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)               &   ! ahmt already * by tmask

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/FLO/flodom.F90

@@ -371 +371 @@
       REAL(wp) ::   zabt, zbct, zcdt, zdat, zabpt, zbcpt, zcdpt, zdapt
       !!---------------------------------------------------------------------
-      !! Statement function
-      REAL(wp) ::   fsline
-      REAL(wp) ::   psax, psay, psbx, psby, psx, psy
-      fsline( psax, psay, psbx, psby, psx, psy ) = psy  * ( psbx - psax )   &
-         &                                       - psx  * ( psby - psay )   &
-         &                                       + psax *   psby - psay * psbx
-      !!---------------------------------------------------------------------
       
       ! 4 semi plane defined by the 4 points and including the T point
@@ -412 +405 @@
    END SUBROUTINE flo_findmesh
 
-
+   FUNCTION fsline( psax, psay, psbx, psby, psx, psy )
+      !! ---------------------------------------------------------------------
+      !!                 ***  Function fsline  ***
+      !!          
+      !! ** Purpose :
+      !! ** Method  : 
+      !!----------------------------------------------------------------------
+      REAL(wp) ::   fsline
+      REAL(wp), INTENT(in) ::   psax, psay, psbx, psby, psx, psy
+      !!---------------------------------------------------------------------
+      fsline = psy  * ( psbx - psax )   &
+         &   - psx  * ( psby - psay )   &
+         &   + psax *   psby - psay * psbx
+      !
+   END FUNCTION fsline
+ 
    FUNCTION flo_dstnce( pla1, phi1, pla2, phi2 )
       !! -------------------------------------------------------------

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/IOM/prtctl.F90

@@ -61 +61 @@
             &                          tab3d_1 = REAL(tab3d_1, 2*wp),                           tab3d_2 =    REAL(tab3d_2, 2*wp), &
             &           mask1 = mask1, mask2 = mask2, &
-            &           clinfo = clinfo, clinfo1 = clinfo1, clinfo2 = clinfo2, clinfo3 = clinfo3 )
+            &           clinfo = clinfo, clinfo1 = clinfo1, clinfo2 = clinfo2, clinfo3 = clinfo3, kdim = kdim )
       ELSEIF( PRESENT(tab2d_1) ) THEN     
          CALL prt_ctl_t(ktab2d_1 = is_tile(tab2d_1), ktab3d_1 = 0, ktab4d_1 = 0, ktab2d_2 = 0, ktab3d_2 = 0,   &
@@ -71 +71 @@
             &                          tab3d_1 =    REAL(tab3d_1, 2*wp),  &
             &           mask1 = mask1,  &
-            &           clinfo = clinfo, clinfo1 = clinfo1, clinfo3 = clinfo3 )
+            &           clinfo = clinfo, clinfo1 = clinfo1, clinfo3 = clinfo3, kdim = kdim )
       ELSEIF( PRESENT(tab4d_1) ) THEN     
          CALL prt_ctl_t(ktab2d_1 = 0, ktab3d_1 = 0, ktab4d_1 = is_tile(tab4d_1), ktab2d_2 = 0, ktab3d_2 = 0,   &
             &                                        tab4d_1 =    REAL(tab4d_1, 2*wp),  &
             &           mask1 = mask1,  &
-            &           clinfo = clinfo, clinfo1 = clinfo1, clinfo3 = clinfo3 )
+            &           clinfo = clinfo, clinfo1 = clinfo1, clinfo3 = clinfo3, kdim = kdim )
       ENDIF
 

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/IOM/restart.F90

@@ -162 +162 @@
       !
       IF( .NOT.ln_diurnal_only ) THEN
-         CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , ssh(:,:        ,Kbb) )     ! before fields
-         CALL iom_rstput( kt, nitrst, numrow, 'ub'     , uu(:,:,:       ,Kbb) )
-         CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vv(:,:,:       ,Kbb) )
-         IF( .NOT.lk_SWE ) THEN
-            CALL iom_rstput( kt, nitrst, numrow, 'tb'  , ts(:,:,:,jp_tem,Kbb) )
-            CALL iom_rstput( kt, nitrst, numrow, 'sb'  , ts(:,:,:,jp_sal,Kbb) )
-         ENDIF
+         CALL iom_rstput( kt, nitrst, numrow, 'sshb', ssh(:,:        ,Kbb) )     ! before fields
+         CALL iom_rstput( kt, nitrst, numrow, 'ub'  , uu(:,:,:       ,Kbb) )
+         CALL iom_rstput( kt, nitrst, numrow, 'vb'  , vv(:,:,:       ,Kbb) )
+         CALL iom_rstput( kt, nitrst, numrow, 'tb'  , ts(:,:,:,jp_tem,Kbb) )
+         CALL iom_rstput( kt, nitrst, numrow, 'sb'  , ts(:,:,:,jp_sal,Kbb) )
          !
 #if ! defined key_RK3
-         CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , ssh(:,:        ,Kmm) )     ! now fields
-         CALL iom_rstput( kt, nitrst, numrow, 'un'     , uu(:,:,:       ,Kmm) )
-         CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vv(:,:,:       ,Kmm) )
-         IF( .NOT.lk_SWE ) THEN
-            CALL iom_rstput( kt, nitrst, numrow, 'tn'  , ts(:,:,:,jp_tem,Kmm) )
-            CALL iom_rstput( kt, nitrst, numrow, 'sn'  , ts(:,:,:,jp_sal,Kmm) )
-            CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop                 )
-         ENDIF
+         CALL iom_rstput( kt, nitrst, numrow, 'sshn', ssh(:,:        ,Kmm) )     ! now fields
+         CALL iom_rstput( kt, nitrst, numrow, 'un'  , uu(:,:,:       ,Kmm) )
+         CALL iom_rstput( kt, nitrst, numrow, 'vn'  , vv(:,:,:       ,Kmm) )
+         CALL iom_rstput( kt, nitrst, numrow, 'tn'  , ts(:,:,:,jp_tem,Kmm) )
+         CALL iom_rstput( kt, nitrst, numrow, 'sn'  , ts(:,:,:,jp_sal,Kmm) )
+         IF( .NOT.lk_SWE )   CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop )
 #endif
       ENDIF
@@ -286 +282 @@
       !                             !*  Read Kbb fields   (NB: in RK3 Kmm = Kbb = Nbb)
       IF(lwp) WRITE(numout,*) '           Kbb u, v and T-S fields read in the restart file'
-      CALL iom_get( numror, jpdom_auto, 'ub'   , uu(:,:,:       ,Kbb), cd_type = 'U', psgn = -1._wp )
-      CALL iom_get( numror, jpdom_auto, 'vb'   , vv(:,:,:       ,Kbb), cd_type = 'V', psgn = -1._wp )
-      IF( .NOT.lk_SWE ) THEN
-         CALL iom_get( numror, jpdom_auto, 'tb', ts(:,:,:,jp_tem,Kbb) )
-         CALL iom_get( numror, jpdom_auto, 'sb', ts(:,:,:,jp_sal,Kbb) )
-      ENDIF
+      CALL iom_get( numror, jpdom_auto, 'ub', uu(:,:,:       ,Kbb), cd_type = 'U', psgn = -1._wp )
+      CALL iom_get( numror, jpdom_auto, 'vb', vv(:,:,:       ,Kbb), cd_type = 'V', psgn = -1._wp )
+      CALL iom_get( numror, jpdom_auto, 'tb', ts(:,:,:,jp_tem,Kbb) )
+      CALL iom_get( numror, jpdom_auto, 'sb', ts(:,:,:,jp_sal,Kbb) )
 #else
       !                             !*  Read Kmm fields   (MLF only)
       IF(lwp) WRITE(numout,*)    '           Kmm u, v and T-S fields read in the restart file'
-      CALL iom_get( numror, jpdom_auto, 'un'   , uu(:,:,:       ,Kmm), cd_type = 'U', psgn = -1._wp )
-      CALL iom_get( numror, jpdom_auto, 'vn'   , vv(:,:,:       ,Kmm), cd_type = 'V', psgn = -1._wp )
-      IF( .NOT.lk_SWE ) THEN
-         CALL iom_get( numror, jpdom_auto, 'tn', ts(:,:,:,jp_tem,Kmm) )
-         CALL iom_get( numror, jpdom_auto, 'sn', ts(:,:,:,jp_sal,Kmm) )
-      ENDIF
+      CALL iom_get( numror, jpdom_auto, 'un', uu(:,:,:       ,Kmm), cd_type = 'U', psgn = -1._wp )
+      CALL iom_get( numror, jpdom_auto, 'vn', vv(:,:,:       ,Kmm), cd_type = 'V', psgn = -1._wp )
+      CALL iom_get( numror, jpdom_auto, 'tn', ts(:,:,:,jp_tem,Kmm) )
+      CALL iom_get( numror, jpdom_auto, 'sn', ts(:,:,:,jp_sal,Kmm) )
       !
       IF( l_1st_euler ) THEN        !*  Euler restart   (MLF only)
@@ -306 +298 @@
          uu(:,:,:  ,Kbb) = uu(:,:,:  ,Kmm)         ! all before fields set to now values
          vv(:,:,:  ,Kbb) = vv(:,:,:  ,Kmm)
-         IF( .NOT.lk_SWE ) ts(:,:,:,:,Kbb) = ts(:,:,:,:,Kmm)
+         ts(:,:,:,:,Kbb) = ts(:,:,:,:,Kmm)
          !
       ELSE                          !* Leap frog restart   (MLF only)
          IF(lwp) WRITE(numout,*) '           Kbb u, v and T-S fields read in the restart file'
-         CALL iom_get( numror, jpdom_auto, 'ub'   , uu(:,:,:       ,Kbb), cd_type = 'U', psgn = -1._wp )
-         CALL iom_get( numror, jpdom_auto, 'vb'   , vv(:,:,:       ,Kbb), cd_type = 'V', psgn = -1._wp )
-         IF( .NOT.lk_SWE ) THEN
-            CALL iom_get( numror, jpdom_auto, 'tb', ts(:,:,:,jp_tem,Kbb) )
-            CALL iom_get( numror, jpdom_auto, 'sb', ts(:,:,:,jp_sal,Kbb) )
-         ENDIF
+         CALL iom_get( numror, jpdom_auto, 'ub', uu(:,:,:       ,Kbb), cd_type = 'U', psgn = -1._wp )
+         CALL iom_get( numror, jpdom_auto, 'vb', vv(:,:,:       ,Kbb), cd_type = 'V', psgn = -1._wp )
+         CALL iom_get( numror, jpdom_auto, 'tb', ts(:,:,:,jp_tem,Kbb) )
+         CALL iom_get( numror, jpdom_auto, 'sb', ts(:,:,:,jp_sal,Kbb) )
       ENDIF
 #endif

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/ISF/isfcav.F90

@@ -170 +170 @@
          ptsc(ji,jj,jp_tem) = zqh(ji,jj) * r1_rho0_rcp
       END_2D
+      CALL lbc_lnk( 'isfmlt', pqfwf, 'T', 1.0_wp)
       !
       ! output fluxes

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/LBC/lbc_lnk_neicoll_generic.h90

@@ -144 +144 @@
       iScnt(:) = PACK( iszall, mask = llsend )                                       ! ok if mask = .false.
       iRcnt(:) = PACK( iszall, mask = llrecv )
-      iSdpl(1) = 0
+      IF( iszS > 0 )   iSdpl(1) = 0
       DO jn = 2,iszS
          iSdpl(jn) = iSdpl(jn-1) + iScnt(jn-1)   ! with _alltoallv: in units of sendtype
       END DO
-      iRdpl(1) = 0
+      IF( iszR > 0 )   iRdpl(1) = 0
       DO jn = 2,iszR
          iRdpl(jn) = iRdpl(jn-1) + iRcnt(jn-1)   ! with _alltoallv: in units of sendtype

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/LBC/lbc_lnk_pt2pt_generic.h90

@@ -123 +123 @@
       !cd     sides:     west  east south north      ;   corners: so-we, so-ea, no-we, no-ea
       isizei(1:4) = (/ ihls, ihls,  ipi,  ipi /)   ;   isizei(5:8) = ihls              ! i- count
-      isizej(1:4) = (/  ipj,  ipj, ihls, ihls /)   ;   isizej(5:8) = ihls              ! j- count
+      isizej(1:4) = (/ Nj_0, Nj_0, ihls, ihls /)   ;   isizej(5:8) = ihls              ! j- count
       ishtSi(1:4) = (/ ip1i, im1i, ip0i, ip0i /)   ;   ishtSi(5:8) = ishtSi( iwewe )   ! i- shift send data
-      ishtSj(1:4) = (/ ip0j, ip0j, ip1j, im1j /)   ;   ishtSj(5:8) = ishtSj( issnn )   ! j- shift send data
+      ishtSj(1:4) = (/ ip1j, ip1j, ip1j, im1j /)   ;   ishtSj(5:8) = ishtSj( issnn )   ! j- shift send data
       ishtRi(1:4) = (/ ip0i, im0i, ip0i, ip0i /)   ;   ishtRi(5:8) = ishtRi( iwewe )   ! i- shift received data location
-      ishtRj(1:4) = (/ ip0j, ip0j, ip0j, im0j /)   ;   ishtRj(5:8) = ishtRj( issnn )   ! j- shift received data location
+      ishtRj(1:4) = (/ ip1j, ip1j, ip0j, im0j /)   ;   ishtRj(5:8) = ishtRj( issnn )   ! j- shift received data location
       ishtPi(1:4) = (/ im1i, ip1i, ip0i, ip0i /)   ;   ishtPi(5:8) = ishtPi( iwewe )   ! i- shift data used for periodicity
-      ishtPj(1:4) = (/ ip0j, ip0j, im1j, ip1j /)   ;   ishtPj(5:8) = ishtPj( issnn )   ! j- shift data used for periodicity
+      ishtPj(1:4) = (/ ip1j, ip1j, im1j, ip1j /)   ;   ishtPj(5:8) = ishtPj( issnn )   ! j- shift data used for periodicity
       !
       ! -------------------------------- !

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/LBC/lbc_nfd_generic.h90

@@ -29 +29 @@
             SELECT CASE ( cd_nat(jf) )
             CASE ( 'T' , 'W' )                         ! T-, W-point
-               DO jl = 1, ipl; DO jk = 1, ipk
+               DO jl = 1, ipl   ;   DO jk = 1, ipk
                   !
                   ! last khls lines (from ipj to ipj-khls+1) : full
@@ -53 +53 @@
                      DO ji = 1, 1                 ! point ipi - khls + 1
                         ii1 = ipi - khls + ji
-                        ii2 =          khls + ji
+                        ii2 =       khls + ji
                         ptab(jf)%pt4d(ii1,ij1,jk,jl) = psgn(jf) * ptab(jf)%pt4d(ii2,ij2,jk,jl)
                      END DO
@@ -84 +84 @@
                END DO; END DO
             CASE ( 'U' )                               ! U-point
-               DO jl = 1, ipl; DO jk = 1, ipk
+               DO jl = 1, ipl   ;   DO jk = 1, ipk
                   !
                   ! last khls lines (from ipj to ipj-khls+1) : full
@@ -103 +103 @@
                      DO ji = 1, khls              ! last khls points
                         ii1 = ipi - khls + ji            ! ends at: ipi - khls + khls = ipi
-                        ii2 = ipi - khls + 1 - ji        ! ends at: ipi - khls + 1 - khls = ipi - 2*khls + 1
+                        ii2 = ipi - khls - ji + 1        ! ends at: ipi - khls + 1 - khls = ipi - 2*khls + 1
                         ptab(jf)%pt4d(ii1,ij1,jk,jl) = psgn(jf) * ptab(jf)%pt4d(ii2,ij2,jk,jl)
                      END DO
@@ -129 +129 @@
                END DO; END DO
             CASE ( 'V' )                               ! V-point
-               DO jl = 1, ipl; DO jk = 1, ipk
+               DO jl = 1, ipl   ;   DO jk = 1, ipk
                   !
                   ! last khls+1 lines (from ipj to ipj-khls) : full
@@ -153 +153 @@
                      DO ji = 1, 1                 ! point ipi - khls + 1
                         ii1 = ipi - khls + ji
-                        ii2 =          khls + ji
+                        ii2 =       khls + ji
                         ptab(jf)%pt4d(ii1,ij1,jk,jl) = psgn(jf) * ptab(jf)%pt4d(ii2,ij2,jk,jl)
                      END DO
@@ -165 +165 @@
                END DO; END DO
             CASE ( 'F' )                               ! F-point
-               DO jl = 1, ipl; DO jk = 1, ipk
+               DO jl = 1, ipl   ;   DO jk = 1, ipk
                   !
                   ! last khls+1 lines (from ipj to ipj-khls) : full
@@ -184 +184 @@
                      DO ji = 1, khls              ! last khls points
                         ii1 = ipi - khls + ji            ! ends at: ipi - khls + khls = ipi
-                        ii2 = ipi - khls + 1 - ji        ! ends at: ipi - khls + 1 - khls = ipi - 2*khls + 1
+                        ii2 = ipi - khls - ji + 1        ! ends at: ipi - khls + 1 - khls = ipi - 2*khls + 1
                         ptab(jf)%pt4d(ii1,ij1,jk,jl) = psgn(jf) * ptab(jf)%pt4d(ii2,ij2,jk,jl)
                      END DO
@@ -198 +198 @@
             SELECT CASE ( cd_nat(jf) )
             CASE ( 'T' , 'W' )                         ! T-, W-point
-               DO jl = 1, ipl; DO jk = 1, ipk
+               DO jl = 1, ipl   ;   DO jk = 1, ipk
                   !
                   ! first: line number ipj-khls : 3 points
@@ -212 +212 @@
                      DO ji = 1, 1                 ! points ipi - khls
                         ii1 = ipi - khls + ji - 1
-                        ii2 =          khls + ji
+                        ii2 =       khls + ji
                         ptab(jf)%pt4d(ii1,ij1,jk,jl) =            ptab(jf)%pt4d(ii2,ij2,jk,jl)   ! Warning: pb with sign...
                      END DO
@@ -240 +240 @@
                      DO ji = 1, khls              ! last khls points
                         ii1 = ipi - khls + ji            ! ends at: ipi - khls + khls = ipi
-                        ii2 = ipi - khls + 1 - ji        ! ends at: ipi - khls + 1 - khls = ipi - 2*khls + 1
+                        ii2 = ipi - khls - ji + 1        ! ends at: ipi - khls + 1 - khls = ipi - 2*khls + 1
                         ptab(jf)%pt4d(ii1,ij1,jk,jl) = psgn(jf) * ptab(jf)%pt4d(ii2,ij2,jk,jl)
                      END DO
@@ -247 +247 @@
                END DO; END DO
             CASE ( 'U' )                               ! U-point
-               DO jl = 1, ipl; DO jk = 1, ipk
+               DO jl = 1, ipl   ;   DO jk = 1, ipk
                   !
                   ! last khls lines (from ipj to ipj-khls+1) : full
@@ -283 +283 @@
                END DO; END DO
             CASE ( 'V' )                               ! V-point
-               DO jl = 1, ipl; DO jk = 1, ipk
+               DO jl = 1, ipl   ;   DO jk = 1, ipk
                   !
                   ! last khls lines (from ipj to ipj-khls+1) : full
@@ -302 +302 @@
                      DO ji = 1, khls            ! last khls points
                         ii1 = ipi - khls + ji          ! ends at: ipi - khls + khls = ipi
-                        ii2 = ipi - khls + 1 - ji      ! ends at: ipi - khls + 1 - khls = ipi - 2*khls + 1
+                        ii2 = ipi - khls - ji + 1      ! ends at: ipi - khls + 1 - khls = ipi - 2*khls + 1
                         ptab(jf)%pt4d(ii1,ij1,jk,jl) = psgn(jf) * ptab(jf)%pt4d(ii2,ij2,jk,jl)
                      END DO
@@ -328 +328 @@
                END DO; END DO
             CASE ( 'F' )                               ! F-point
-               DO jl = 1, ipl; DO jk = 1, ipk
+               DO jl = 1, ipl   ;   DO jk = 1, ipk
                   !
                   ! last khls lines (from ipj to ipj-khls+1) : full

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/LBC/lbcnfd.F90

@@ -10 +10 @@
 
    !!----------------------------------------------------------------------
-   !!   lbc_nfd       : generic interface for lbc_nfd_3d and lbc_nfd_2d routines
-   !!   lbc_nfd_3d    : lateral boundary condition: North fold treatment for a 3D arrays   (lbc_nfd)
-   !!   lbc_nfd_2d    : lateral boundary condition: North fold treatment for a 2D arrays   (lbc_nfd)
-   !!   lbc_nfd_nogather       : generic interface for lbc_nfd_nogather_3d and 
-   !!                            lbc_nfd_nogather_2d routines (designed for use
-   !!                            with ln_nnogather to avoid global width arrays
-   !!                            mpi all gather operations)
+   !!   lbc_nfd       : generic interface for lbc_nfd_sp and lbc_nfd_dp routines that is doing the north fold in a non-mpi case 
+   !!   mpp_nfd       : generic interface for mpp_nfd_sp and mpp_nfd_dp routines that will use lbc_nfd directly or indirectly
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain 
@@ -36 +31 @@
       MODULE PROCEDURE   mpp_nfd_sp, mpp_nfd_dp
    END INTERFACE
-
-   INTERFACE lbc_nfd_nogather   ! called by mpp_nfd
-      MODULE PROCEDURE   lbc_nfd_nogather_sp, lbc_nfd_nogather_dp
-   END INTERFACE
    
    PUBLIC   mpp_nfd            ! mpi north fold conditions
    PUBLIC   lbc_nfd            ! north fold conditions
-   PUBLIC   lbc_nfd_nogather   ! north fold conditions (no allgather case)
 
-   INTEGER, PUBLIC, PARAMETER            ::   jpmaxngh = 3               !:
-   INTEGER, PUBLIC                       ::   nsndto                     !:
-   INTEGER, PUBLIC, DIMENSION (jpmaxngh) ::   isendto                    !: processes to which communicate
-   INTEGER, PUBLIC                       ::   ijpj
-
+   INTEGER, PUBLIC                               :: nfd_nbnei
+   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION (:  ) :: nfd_rknei
+   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION (:,:) :: nfd_rksnd
+   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION (:,:) :: nfd_jisnd
+   
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -59 +49 @@
    !!----------------------------------------------------------------------
    !!                   ***  routine lbc_nfd_[sd]p  ***
-   !!               ***  routine lbc_nfd_nogather_[sd]p  ***
    !!                   ***  routine lbc_nfd_ext_[sd]p  ***
    !!----------------------------------------------------------------------
@@ -75 +64 @@
 #define PRECISION sp
 #  include "lbc_nfd_generic.h90"
-#  include "lbc_nfd_nogather_generic.h90"
 #  include "lbc_nfd_ext_generic.h90"
 #undef PRECISION
@@ -83 +71 @@
 #define PRECISION dp
 #  include "lbc_nfd_generic.h90"
-#  include "lbc_nfd_nogather_generic.h90"
 #  include "lbc_nfd_ext_generic.h90"
 #undef PRECISION

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/LBC/lib_mpp.F90

@@ -49 +49 @@
    !!   mppsync       :
    !!   mppstop       :
-   !!   mpp_ini_north : initialisation of north fold
+   !!   mpp_ini_northgather : initialisation of north fold with gathering of the communications
    !!   mpp_lbc_north_icb : alternative to mpp_nfd for extra outer halo with icebergs
    !!   mpp_bcast_nml : broadcast/receive namelist character buffer from reading process to all others
@@ -64 +64 @@
    PUBLIC   ctl_stop, ctl_warn, ctl_opn, ctl_nam, load_nml
    PUBLIC   mpp_start, mppstop, mppsync, mpp_comm_free
-   PUBLIC   mpp_ini_north
+   PUBLIC   mpp_ini_northgather
    PUBLIC   mpp_min, mpp_max, mpp_sum, mpp_minloc, mpp_maxloc
    PUBLIC   mpp_delay_max, mpp_delay_sum, mpp_delay_rcv
@@ -1152 +1152 @@
 
 
-   SUBROUTINE mpp_ini_north
-      !!----------------------------------------------------------------------
-      !!               ***  routine mpp_ini_north  ***
+   SUBROUTINE mpp_ini_northgather
+      !!----------------------------------------------------------------------
+      !!               ***  routine mpp_ini_northgather  ***
       !!
       !! ** Purpose :   Initialize special communicator for north folding
@@ -1210 +1210 @@
       !
 #endif
-   END SUBROUTINE mpp_ini_north
+   END SUBROUTINE mpp_ini_northgather
 
 

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/LBC/mpp_nfd_generic.h90

@@ -10 +10 @@
       !
       LOGICAL  ::   ll_add_line
-      INTEGER  ::   ji,  jj,  jk,  jl, jh, jf, jr   ! dummy loop indices
+      INTEGER  ::   ji,  jj,  jk,  jl, jf, jr, jg, jn   ! dummy loop indices
       INTEGER  ::   ipi, ipj, ipj2, ipk, ipl, ipf   ! dimension of the input array
-      INTEGER  ::   imigr, iihom, ijhom             ! local integers
       INTEGER  ::   ierr, ibuffsize, iis0, iie0, impp
-      INTEGER  ::   ii1, ii2, ij1, ij2
-      INTEGER  ::   ipimax, i0max
+      INTEGER  ::   ii1, ii2, ij1, ij2, iis, iie, iib, iig, iin
+      INTEGER  ::   i0max
       INTEGER  ::   ij, iproc, ipni, ijnr
-      INTEGER, DIMENSION (jpmaxngh)       ::   ml_req_nf   ! for mpi_isend when avoiding mpi_allgather
-      INTEGER                             ::   ml_err      ! for mpi_isend when avoiding mpi_allgather
-      !                                                    ! Workspace for message transfers avoiding mpi_allgather
-      INTEGER                             ::   ipj_b       ! sum of lines for all multi fields
-      INTEGER                             ::   i012        ! 0, 1 or 2
-      INTEGER , DIMENSION(:,:)        , ALLOCATABLE ::   jj_s  ! position of sent lines
-      INTEGER , DIMENSION(:,:)        , ALLOCATABLE ::   jj_b  ! position of buffer lines
-      INTEGER , DIMENSION(:)          , ALLOCATABLE ::   ipj_s ! number of sent lines
-      REAL(PRECISION), DIMENSION(:,:,:,:)    , ALLOCATABLE ::   ztabb, ztabr, ztabw  ! buffer, receive and work arrays
+      INTEGER, DIMENSION (:), ALLOCATABLE ::   ireq_s, ireq_r   ! for mpi_isend when avoiding mpi_allgather
+      INTEGER                             ::   ipjtot           ! sum of lines for all multi fields
+      INTEGER                             ::   i012             ! 0, 1 or 2
+      INTEGER , DIMENSION(:,:)        , ALLOCATABLE ::   ijsnd  ! j-position of sent lines for each field
+      INTEGER , DIMENSION(:,:)        , ALLOCATABLE ::   ijbuf  ! j-position of send buffer lines for each field
+      INTEGER , DIMENSION(:,:)        , ALLOCATABLE ::   ijrcv  ! j-position of recv buffer lines for each field
+      INTEGER , DIMENSION(:,:)        , ALLOCATABLE ::   ii1st, iiend
+      INTEGER , DIMENSION(:)          , ALLOCATABLE ::   ipjfld ! number of sent lines for each field
+      REAL(PRECISION), DIMENSION(:,:,:,:)    , ALLOCATABLE ::   zbufs  ! buffer, receive and work arrays
+      REAL(PRECISION), DIMENSION(:,:,:,:,:)  , ALLOCATABLE ::   zbufr  ! buffer, receive and work arrays
       REAL(PRECISION), DIMENSION(:,:,:,:,:)  , ALLOCATABLE ::   znorthloc
       REAL(PRECISION), DIMENSION(:,:,:,:,:,:), ALLOCATABLE ::   znorthglo
@@ -62 +62 @@
          ll_add_line = l_full_nf_update .OR. ( ncom_stp <= nit000+1 .AND. .NOT. ln_rstart )
          
-         ALLOCATE(ipj_s(ipf))                ! how many lines do we exchange?
+         ALLOCATE(ipjfld(ipf))                 ! how many lines do we exchange for each field?
          IF( ll_add_line ) THEN
-            DO jf = 1, ipf                      ! Loop over the number of arrays to be processed
-               ipj_s(jf) = khls + COUNT( (/ c_NFtype == 'T' .OR. cd_nat(jf) == 'V' .OR. cd_nat(jf) == 'F' /) ) 
+            DO jf = 1, ipf                     ! Loop over the number of arrays to be processed
+               ipjfld(jf) = khls + COUNT( (/ c_NFtype == 'T' .OR. cd_nat(jf) == 'V' .OR. cd_nat(jf) == 'F' /) )
             END DO
          ELSE
-            ipj_s(:) = khls
+            ipjfld(:) = khls
          ENDIF
          
-         ipj   = MAXVAL(ipj_s(:))            ! Max 2nd dimension of message transfers
-         ipj_b = SUM(   ipj_s(:))            ! Total number of lines to be exchanged
-         ALLOCATE( jj_s(ipj, ipf), jj_b(ipj, ipf) )
+         ipj    = MAXVAL(ipjfld(:))            ! Max 2nd dimension of message transfers
+         ipjtot = SUM(   ipjfld(:))            ! Total number of lines to be exchanged
 
          ! Index of modifying lines in input
+         ALLOCATE( ijsnd(ipj, ipf), ijbuf(ipj, ipf), ijrcv(ipj, ipf), ii1st(ipj, ipf), iiend(ipj, ipf) )
+
          ij1 = 0
-         DO jf = 1, ipf                      ! Loop over the number of arrays to be processed
-            !
+         DO jf = 1, ipf                        ! Loop over the number of arrays to be processed
+            !
+            DO jj = 1, khls   ! first khls lines (starting from top) must be fully defined
+               ii1st(jj, jf) = 1
+               iiend(jj, jf) = jpi
+            END DO
+            !
+            ! what do we do with line khls+1 (starting from top)
             IF( c_NFtype == 'T' ) THEN          ! *  North fold  T-point pivot
                SELECT CASE ( cd_nat(jf) )
-               CASE ( 'T', 'W', 'U' )   ;   i012 = 1   ! T-, U-, W-point
-               CASE ( 'V', 'F'      )   ;   i012 = 2   ! V-, F-point
+               CASE ('T','W')   ;   i012 = 1   ;   ii1st(khls+1, jf) = mi0(jpiglo/2+2)   ;   iiend(khls+1, jf) = mi1(jpiglo-khls)
+               CASE ('U'    )   ;   i012 = 1   ;   ii1st(khls+1, jf) = mi0(jpiglo/2+1)   ;   iiend(khls+1, jf) = mi1(jpiglo-khls)
+               CASE ('V'    )   ;   i012 = 2   ;   ii1st(khls+1, jf) = 1                 ;   iiend(khls+1, jf) = jpi
+               CASE ('F'    )   ;   i012 = 2   ;   ii1st(khls+1, jf) = 1                 ;   iiend(khls+1, jf) = jpi
                END SELECT
             ENDIF
             IF( c_NFtype == 'F' ) THEN          ! *  North fold  F-point pivot
                SELECT CASE ( cd_nat(jf) )
-               CASE ( 'T', 'W', 'U' )   ;   i012 = 0   ! T-, U-, W-point
-               CASE ( 'V', 'F'      )   ;   i012 = 1   ! V-, F-point
+               CASE ('T','W')   ;   i012 = 0   ! we don't touch line khls+1
+               CASE ('U'    )   ;   i012 = 0   ! we don't touch line khls+1
+               CASE ('V'    )   ;   i012 = 1   ;   ii1st(khls+1, jf) = mi0(jpiglo/2+1)   ;   iiend(khls+1, jf) = mi1(jpiglo-khls  )
+               CASE ('F'    )   ;   i012 = 1   ;   ii1st(khls+1, jf) = mi0(jpiglo/2+1)   ;   iiend(khls+1, jf) = mi1(jpiglo-khls-1)
                END SELECT
             ENDIF
-               !
-            DO jj = 1, ipj_s(jf)
+            !
+            DO jj = 1, ipjfld(jf)
                ij1 = ij1 + 1
-               jj_b(jj,jf) = ij1
-               jj_s(jj,jf) = jpj - 2*khls + jj - i012
+               ijsnd(jj,jf) = jpj - 2*khls + jj - i012   ! sent lines (from bottom of sent lines)
+               ijbuf(jj,jf) = ij1                        ! gather all lines in the snd/rcv buffers
+               ijrcv(jj,jf) = jpj - jj + 1               ! recv lines (from the top -> reverse order for jj)
             END DO
             !
          END DO
          !
-         ALLOCATE( ztabb(jpimax,ipj_b,ipk,ipl) )   ! store all the data to be sent in a buffer array
-         ibuffsize = jpimax * ipj_b * ipk * ipl
-         !
+         i0max = jpimax - 2 * khls                                    ! we are not sending the halos
+         ALLOCATE( zbufs(i0max,ipjtot,ipk,ipl), ireq_s(nfd_nbnei) )   ! store all the data to be sent in a buffer array
+         ibuffsize = i0max * ipjtot * ipk * ipl
+         !
+         ! fill the send buffer with all the lines
          DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk
-            DO jj = 1, ipj_s(jf)
-               ij1 = jj_b(jj,jf)
-               ij2 = jj_s(jj,jf)
-               DO ji = 1, jpi
-                  ztabb(ji,ij1,jk,jl) = ptab(jf)%pt4d(ji,ij2,jk,jl)
-               END DO
-               DO ji = jpi+1, jpimax
-                  ztabb(ji,ij1,jk,jl) = HUGE(0._/**/PRECISION)   ! avoid sending uninitialized values (make sure we don't use it)
+            DO jj = 1, ipjfld(jf)
+               ij1 = ijbuf(jj,jf)
+               ij2 = ijsnd(jj,jf)
+               DO ji = Nis0, Nie0       ! should not use any other value
+                  iib = ji - Nis0 + 1
+                  zbufs(iib,ij1,jk,jl) = ptab(jf)%pt4d(ji,ij2,jk,jl)
+               END DO
+               DO ji = Ni_0+1, i0max    ! avoid sending uninitialized values (make sure we don't use it)
+                  zbufs(ji,ij1,jk,jl) = HUGE(0._/**/PRECISION)   ! make sure we don't use it...
                END DO
             END DO
@@ -119 +134 @@
          IF( ln_timing ) CALL tic_tac(.TRUE.)
          !
-         ! send the data as soon as possible
-         DO jr = 1, nsndto
-            iproc = nfproc(isendto(jr))
+         ! send the same buffer data to all neighbourgs as soon as possible
+         DO jn = 1, nfd_nbnei
+            iproc = nfd_rknei(jn)
             IF( iproc /= narea-1 .AND. iproc /= -1 ) THEN
 #if ! defined key_mpi_off
-               CALL MPI_ISEND( ztabb, ibuffsize, MPI_TYPE, iproc, 5, mpi_comm_oce, ml_req_nf(jr), ierr )
+               CALL MPI_Isend( zbufs, ibuffsize, MPI_TYPE, iproc, 5, mpi_comm_oce, ireq_s(jn), ierr )
 #endif
+            ELSE
+               ireq_s(jn) = MPI_REQUEST_NULL
             ENDIF
          END DO
          !
-         ipimax = jpimax * jpmaxngh
-         ALLOCATE( ztabw(jpimax,ipj_b,ipk,ipl), ztabr(ipimax,ipj_b,ipk,ipl) ) 
-         !
-         DO jr = 1, nsndto
-            !
-            ipni  = isendto(jr)
-            iproc = nfproc(ipni)
-            ipi   = nfjpi (ipni)
-            !
-            IF( ipni ==   1  ) THEN   ;   iis0 =   1          ! domain  left side: as e-w comm already done -> from 1st column
-            ELSE                      ;   iis0 =   1 + khls   ! default: -> from inner domain 
-            ENDIF
-            IF( ipni == jpni ) THEN   ;   iie0 = ipi          ! domain right side: as e-w comm already done -> until last column
-            ELSE                      ;   iie0 = ipi - khls   ! default: -> until inner domain 
-            ENDIF
-            impp = nfimpp(ipni) - nfimpp(isendto(1))
+         ALLOCATE( zbufr(i0max,ipjtot,ipk,ipl,nfd_nbnei), ireq_r(nfd_nbnei) ) 
+         !
+         DO jn = 1, nfd_nbnei
+            !
+            iproc = nfd_rknei(jn)
             !
             IF(           iproc == -1 ) THEN   ! No neighbour (land proc that was suppressed)
                !
+               ireq_r(jn) = MPI_REQUEST_NULL                ! no message to be received
+               zbufr(:,:,:,:,jn) = HUGE(0._/**/PRECISION)   ! default: define it and make sure we don't use it...
                SELECT CASE ( kfillmode )
-               CASE ( jpfillnothing )               ! no filling 
-               CASE ( jpfillcopy    )               ! filling with inner domain values
+               CASE ( jpfillnothing )                       ! no filling 
+               CASE ( jpfillcopy    )                       ! filling with inner domain values
                   DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk
-                     DO jj = 1, ipj_s(jf)
-                        ij1 = jj_b(jj,jf)
-                        ij2 = jj_s(jj,jf)
-                        DO ji = iis0, iie0
-                           ztabr(impp+ji,ij1,jk,jl) = ptab(jf)%pt4d(Nis0,ij2,jk,jl)   ! chose to take the 1st iner domain point
-                        END DO
+                     DO jj = 1, ipjfld(jf)
+                        ij1 = ijbuf(jj,jf)
+                        ij2 = ijsnd(jj,jf)                                      ! we will use only the first value, see init_nfdcom
+                        zbufr(1,ij1,jk,jl,jn) = ptab(jf)%pt4d(Nis0,ij2,jk,jl)   ! chose to take the 1st inner domain point
                      END DO
                   END DO   ;   END DO   ;   END DO
-               CASE ( jpfillcst     )               ! filling with constant value
-                  DO jl = 1, ipl   ;   DO jk = 1, ipk
-                     DO jj = 1, ipj_b
-                        DO ji = iis0, iie0
-                           ztabr(impp+ji,jj,jk,jl) = pfillval
-                        END DO
-                     END DO
-                  END DO   ;   END DO
+               CASE ( jpfillcst     )                       ! filling with constant value
+                  zbufr(1,:,:,:,jn) = pfillval              ! we will use only the first value, see init_nfdcom
                END SELECT
                !
             ELSE IF( iproc == narea-1 ) THEN   ! get data from myself!
                !
+               ireq_r(jn) = MPI_REQUEST_NULL                ! no message to be received
                DO jf = 1, ipf   ;   DO jl = 1, ipl  ;   DO jk = 1, ipk
-                  DO jj = 1, ipj_s(jf)
-                     ij1 = jj_b(jj,jf)
-                     ij2 = jj_s(jj,jf)
-                     DO ji = iis0, iie0
-                        ztabr(impp+ji,ij1,jk,jl) = ptab(jf)%pt4d(ji,ij2,jk,jl)
+                  DO jj = 1, ipjfld(jf)
+                     ij1 = ijbuf(jj,jf)
+                     ij2 = ijsnd(jj,jf)
+                     DO ji = Nis0, Nie0                     ! should not use any other value
+                        iib = ji - Nis0 + 1
+                        zbufr(iib,ij1,jk,jl,jn) = ptab(jf)%pt4d(ji,ij2,jk,jl)
                      END DO
                   END DO
@@ -183 +185 @@
                !
             ELSE                               ! get data from a neighbour trough communication
-               !  
 #if ! defined key_mpi_off
-               CALL MPI_RECV( ztabw, ibuffsize, MPI_TYPE, iproc, 5, mpi_comm_oce, MPI_STATUS_IGNORE, ierr )
+               CALL MPI_Irecv( zbufr(:,:,:,:,jn), ibuffsize, MPI_TYPE, iproc, 5, mpi_comm_oce, ireq_r(jn), ierr )
 #endif
-               DO jl = 1, ipl   ;   DO jk = 1, ipk
-                  DO jj = 1, ipj_b
-                     DO ji = iis0, iie0
-                        ztabr(impp+ji,jj,jk,jl) = ztabw(ji,jj,jk,jl)
-                     END DO
-                  END DO
-               END DO   ;   END DO
-               
-            ENDIF
-            !
-         END DO   ! nsndto
+            ENDIF
+            !
+         END DO   ! nfd_nbnei
+         !
+         CALL mpi_waitall(nfd_nbnei, ireq_r, MPI_STATUSES_IGNORE, ierr)   ! wait for all Irecv
          !
          IF( ln_timing ) CALL tic_tac(.FALSE.)
@@ -204 +199 @@
          !
          DO jf = 1, ipf
-            ij1 = jj_b(       1 ,jf)
-            ij2 = jj_b(ipj_s(jf),jf)
-            CALL lbc_nfd_nogather( ptab(jf)%pt4d(:,:,:,:), ztabr(:,ij1:ij2,:,:), cd_nat(jf), psgn(jf), khls )
-         END DO
-         !
-         DEALLOCATE( ztabr, ztabw, jj_s, jj_b, ipj_s )
-         !
-         DO jr = 1,nsndto
-            iproc = nfproc(isendto(jr))
-            IF( iproc /= narea-1 .AND. iproc /= -1 ) THEN
-               CALL mpi_wait( ml_req_nf(jr), MPI_STATUS_IGNORE, ml_err )   ! put the wait at the very end just before the deallocate
-            ENDIF
-         END DO
-         DEALLOCATE( ztabb )
+            !
+            SELECT CASE ( cd_nat(jf) )     ! which grid number?
+            CASE ('T','W')   ;   iig = 1   ! T-, W-point
+            CASE ('U')       ;   iig = 2   ! U-point
+            CASE ('V')       ;   iig = 3   ! V-point
+            CASE ('F')       ;   iig = 4   ! F-point
+            END SELECT
+            !
+            DO jl = 1, ipl   ;   DO jk = 1, ipk
+               !
+               ! if T point with F-point pivot : must be done first
+               !    --> specific correction of 3 points near the 2 pivots (to be clean, usually masked -> so useless) 
+               IF( c_NFtype == 'F' .AND. iig == 1 ) THEN
+                  ij1 = jpj - khls     ! j-index in the receiving array
+                  ij2 = 1              ! only 1 line in the buffer
+                  DO ji = mi0(khls), mi1(khls)
+                     iib = nfd_jisnd(mi0(       khls),iig)   ! i-index in the buffer
+                     iin = nfd_rksnd(mi0(       khls),iig)   ! neigbhour-index in the buffer
+                     IF( nfd_rknei(iin) == -1 .AND. kfillmode == jpfillnothing )   CYCLE
+                     ptab(jf)%pt4d(ji,ij1,jk,jl) = zbufr(iib,ij2,jk,jl,iin)   ! no psgn(jf)
+                  END DO
+                  DO ji = mi0(jpiglo/2+1), mi1(jpiglo/2+1)
+                     iib = nfd_jisnd(mi0( jpiglo/2+1),iig)   ! i-index in the buffer
+                     iin = nfd_rksnd(mi0( jpiglo/2+1),iig)   ! neigbhour-index in the buffer
+                     IF( nfd_rknei(iin) == -1 .AND. kfillmode == jpfillnothing )   CYCLE
+                     ptab(jf)%pt4d(ji,ij1,jk,jl) = zbufr(iib,ij2,jk,jl,iin)   ! no psgn(jf)
+                  END DO
+                  DO ji = mi0(jpiglo-khls), mi1(jpiglo-khls)
+                     iib = nfd_jisnd(mi0(jpiglo-khls),iig)   ! i-index in the buffer
+                     iin = nfd_rksnd(mi0(jpiglo-khls),iig)   ! neigbhour-index in the buffer
+                     IF( nfd_rknei(iin) == -1 .AND. kfillmode == jpfillnothing )   CYCLE
+                     ptab(jf)%pt4d(ji,ij1,jk,jl) = zbufr(iib,ij2,jk,jl,iin)   ! no psgn(jf)
+                  END DO
+               ENDIF
+               !
+               ! Apply the North pole folding.
+               DO jj = 1, ipjfld(jf)   ! for all lines to be exchanged for this field
+                  ij1 = ijrcv(jj,jf)   ! j-index in the receiving array
+                  ij2 = ijbuf(jj,jf)   ! j-index in the buffer
+                  iis = ii1st(jj,jf)   ! stating i-index in the receiving array
+                  iie = iiend(jj,jf)   !  ending i-index in the receiving array
+                  DO ji = iis, iie 
+                     iib = nfd_jisnd(ji,iig)   ! i-index in the buffer
+                     iin = nfd_rksnd(ji,iig)   ! neigbhour-index in the buffer
+                     IF( nfd_rknei(iin) == -1 .AND. kfillmode == jpfillnothing )   CYCLE
+                     ptab(jf)%pt4d(ji,ij1,jk,jl) = psgn(jf) * zbufr(iib,ij2,jk,jl,iin)
+                  END DO
+               END DO
+               !
+               ! re-apply periodocity when we modified the eastern side of the inner domain (and not the full line)
+               IF( c_NFtype == 'T' ) THEN          ! *  North fold  T-point pivot
+                  IF(     iig <= 2 ) THEN   ;   iis = mi0(1)   ;   iie = mi1(khls)   ! 'T','W','U': update west halo
+                  ELSE                      ;   iis = 1        ;   iie = 0           ! 'V','F'    : full line already exchanged
+                  ENDIF
+               ENDIF
+               IF( c_NFtype == 'F' ) THEN          ! *  North fold  F-point pivot
+                  IF(     iig <= 2 ) THEN   ;   iis = 1        ;   iie = 0           ! 'T','W','U': nothing to do
+                  ELSEIF( iig == 3 ) THEN   ;   iis = mi0(1)   ;   iie = mi1(khls)   ! 'V'        : update west halo
+                  ELSEIF( khls > 1 ) THEN   ;   iis = mi0(1)   ;   iie = mi1(khls-1) ! 'F' and khls > 1
+                  ELSE                      ;   iis = 1        ;   iie = 0           ! 'F' and khls == 1 : nothing to do
+                  ENDIF
+               ENDIF
+               jj  = ipjfld(jf)     ! only for the last line of this field
+               ij1 = ijrcv(jj,jf)   ! j-index in the receiving array
+               ij2 = ijbuf(jj,jf)   ! j-index in the buffer
+               DO ji = iis, iie
+                  iib = nfd_jisnd(ji,iig)   ! i-index in the buffer
+                  iin = nfd_rksnd(ji,iig)   ! neigbhour-index in the buffer
+                  IF( nfd_rknei(iin) == -1 .AND. kfillmode == jpfillnothing )   CYCLE
+                  ptab(jf)%pt4d(ji,ij1,jk,jl) = psgn(jf) * zbufr(iib,ij2,jk,jl,iin)
+               END DO
+               !               
+            END DO   ;   END DO   ! ipl   ; ipk
+            !               
+         END DO   ! ipf
+       
+         !
+         DEALLOCATE( zbufr, ireq_r, ijsnd, ijbuf, ijrcv, ii1st, iiend, ipjfld )
+         !
+         CALL mpi_waitall(nfd_nbnei, ireq_s, MPI_STATUSES_IGNORE, ierr)   ! wait for all Isend
+         !
+         DEALLOCATE( zbufs, ireq_s )
          !
       ELSE                             !==  allgather exchanges  ==!
@@ -265 +328 @@
               !
                SELECT CASE ( kfillmode )
-               CASE ( jpfillnothing )               ! no filling 
+               CASE ( jpfillnothing )               ! no filling
+                  CALL ctl_stop( 'STOP', 'mpp_nfd_generic : cannot use jpfillnothing with ln_nnogather = F')
                CASE ( jpfillcopy    )               ! filling with inner domain values
                   DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk
@@ -329 +393 @@
          DEALLOCATE( ztabglo )
          !
-      ENDIF   ! l_north_nogather
+      ENDIF   ! ln_nnogather
       !
    END SUBROUTINE mpp_nfd_/**/PRECISION

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/LBC/mppini.F90

@@ -23 +23 @@
    USE bdy_oce        ! open BounDarY
    !
-   USE lbcnfd  , ONLY : isendto, nsndto ! Setup of north fold exchanges
+   USE lbcnfd         ! Setup of north fold exchanges
    USE lib_mpp        ! distribued memory computing library
    USE iom            ! nemo I/O library
@@ -168 +168 @@
       IF(lwp) THEN
             WRITE(numout,*) '   Namelist nammpp'
-         IF( jpni < 1 .OR. jpnj < 1  ) THEN
+         IF( jpni < 1 .OR. jpnj < 1 ) THEN
             WRITE(numout,*) '      jpni and jpnj will be calculated automatically'
          ELSE
@@ -331 +331 @@
       njmpp = ijmppt(ii,ij)
       !
-      CALL init_doloop                          ! set start/end indices of do-loop, depending on the halo width value (nn_hls)
+      CALL init_doloop    ! set start/end indices of do-loop, depending on the halo width value (nn_hls)
+      CALL init_locglo    ! define now functions needed to convert indices from/to global to/from local domains
       !
       IF(lwp) THEN
@@ -503 +504 @@
          WRITE(numout,*) '      mpi nei no-we = ', mpinei(jpnw)  , '   mpi nei no-ea = ', mpinei(jpne)
       ENDIF
-      !                          ! Prepare mpp north fold
-      !
-      llmpiNFold =          jpni  > 1 .AND. l_NFold   ! is the North fold done with an MPI communication?
-      l_IdoNFold = ijn(narea) == jpnj .AND. l_NFold   ! is this process doing North fold?
-      !
-      IF( llmpiNFold ) THEN
-         CALL mpp_ini_north
-         IF (lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) '   ==>>>   North fold boundary prepared for jpni >1'
-         ENDIF
-         IF (llwrtlay) THEN      ! additional prints in layout.dat
-            WRITE(inum,*)
-            WRITE(inum,*)
-            WRITE(inum,*) 'Number of subdomains located along the north fold : ', ndim_rank_north
-            WRITE(inum,*) 'Rank of the subdomains located along the north fold : ', ndim_rank_north
-            DO jp = 1, ndim_rank_north, 5
-               WRITE(inum,*) nrank_north( jp:MINVAL( (/jp+4,ndim_rank_north/) ) )
-            END DO
-         ENDIF
-         IF ( l_IdoNFold .AND. ln_nnogather ) THEN
-            CALL init_nfdcom     ! northfold neighbour lists
-            IF (llwrtlay) THEN
-               WRITE(inum,*)
-               WRITE(inum,*) 'north fold exchanges with explicit point-to-point messaging :'
-               WRITE(inum,*) '   nsndto  : ', nsndto
-               WRITE(inum,*) '   isendto : ', isendto(1:nsndto)
-            ENDIF
-         ENDIF
-      ENDIF
       !
       CALL mpp_ini_nc(nn_hls)    ! Initialize communicator for neighbourhood collective communications
@@ -539 +510 @@
          mpi_nc_com8(jh) = mpi_nc_com8(nn_hls)
       END DO
-      !
-      IF( jpnij > 1) CALL init_excl_landpt      ! exclude exchanges which contain only land points
-      !
-      ! Save processor layout changes in ascii file
+      !                          ! Exclude exchanges which contain only land points
+      !
+      IF( jpnij > 1 ) CALL init_excl_landpt
+      !
+      !                          ! Prepare mpp north fold
+      !
+      llmpiNFold =          jpni  > 1 .AND. l_NFold           ! is the North fold done with an MPI communication?
+      l_IdoNFold = ijn(narea) == jpnj .AND. l_NFold           ! is this process doing North fold?
+      !
+      IF( llmpiNFold )   CALL init_nfdcom( llwrtlay, inum )   ! init northfold communication, must be done after init_excl_landpt
+      !
+      !                          !  Save processor layout changes in ascii file
+      !
       DO jh = 1, n_hlsmax    ! different halo size
          DO ji = 1, 8
@@ -732 +712 @@
       isziref = jpiglo*jpjglo+1   ! define a value that is larger than the largest possible
       iszjref = jpiglo*jpjglo+1
+      !
+      ! WARNING, see also init_excl_landpt: minimum subdomain size defined here according to nn_hls (and not n_hlsmax)
+      ! --> If, one day, we want to use local halos largers than nn_hls, we must replace nn_hls by n_hlsmax
       !
       iszimin = 3*nn_hls          ! minimum size of the MPI subdomain so halos are always adressing neighbor inner domain
@@ -1165 +1148 @@
       INTEGER ::   iiwe, iiea, iist, iisz 
       INTEGER ::   ijso, ijno, ijst, ijsz 
-      LOGICAL ::   llsave
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zmsk
       LOGICAL , DIMENSION(Ni_0,Nj_0,1)      ::   lloce
@@ -1174 +1156 @@
       !
       ! Here we look only at communications excluding the NP folding.
-      ! As lbcnfd not validated if halo size /= nn_hls, we switch if off temporary...
-      llsave = l_IdoNFold
+      !   --> we switch off lbcnfd at this stage (init_nfdcom called after init_excl_landpt)...
       l_IdoNFold = .FALSE.
       !
-      DO jh = 1, n_hlsmax    ! different halo size
+      ! WARNING, see also bestpartition: minimum subdomain size defined in bestpartition according to nn_hls.
+      ! If, one day, we want to use local halos largers than nn_hls, we must replace nn_hls by n_hlsmax in bestpartition
+      !
+      DO jh = 1, MIN(nn_hls, n_hlsmax)   ! different halo size
          !
          ipi = Ni_0 + 2*jh   ! local domain size
@@ -1186 +1170 @@
          zmsk(jh+1:jh+Ni_0,jh+1:jh+Nj_0) = REAL(COUNT(lloce, dim = 3), wp)   ! define inner domain -> need REAL to use lbclnk
          CALL lbc_lnk('mppini', zmsk, 'T', 1._wp, khls = jh)                 ! fill halos
+         ! Beware, coastal F points can be used in the code -> we may need communications for these points F points even if tmask = 0
+         ! -> the mask we must use here is equal to 1 as soon as one of the 4 neighbours is oce (sum of the mask, not multiplication)
+         zmsk(jh+1:jh+Ni_0,jh+1:jh+Nj_0) = zmsk(jh+1:jh+Ni_0,jh+1  :jh+Nj_0  ) + zmsk(jh+1+1:jh+Ni_0+1,jh+1  :jh+Nj_0  )   &
+            &                            + zmsk(jh+1:jh+Ni_0,jh+1+1:jh+Nj_0+1) + zmsk(jh+1+1:jh+Ni_0+1,jh+1+1:jh+Nj_0+1)
+         CALL lbc_lnk('mppini', zmsk, 'T', 1._wp, khls = jh)                 ! fill halos again!
          !        
-         iiwe = jh   ;   iiea = Ni_0   ! bottom-left corfer - 1 of the sent data
+         iiwe = jh   ;   iiea = Ni_0   ! bottom-left corner - 1 of the sent data
          ijso = jh   ;   ijno = Nj_0
          IF( nn_comm == 1 ) THEN 
             iist =  0   ;   iisz = ipi
-            ijst =  0   ;   ijsz = ipj
+            ijst = jh   ;   ijsz = Nj_0
          ELSE
             iist = jh   ;   iisz = Ni_0
@@ -1207 +1196 @@
          IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijno+1:ijno+jh  ) )) == 0 )   mpiSnei(jh,jpne) = -1
          !
-         iiwe = iiwe-jh   ;   iiea = iiea+jh   ! bottom-left corfer - 1 of the received data
+         iiwe = iiwe-jh   ;   iiea = iiea+jh   ! bottom-left corner - 1 of the received data
          ijso = ijso-jh   ;   ijno = ijno+jh
          ! do not send if we send only land points
@@ -1233 +1222 @@
          !
       END DO
-      l_IdoNFold = llsave
 
    END SUBROUTINE init_excl_landpt
@@ -1278 +1266 @@
 
 
-   SUBROUTINE init_nfdcom
+   SUBROUTINE init_nfdcom( ldwrtlay, knum )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE  init_nfdcom  ***
@@ -1288 +1276 @@
       !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)
       !!    2.0  ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
-      !!----------------------------------------------------------------------
-      INTEGER  ::   sxM, dxM, sxT, dxT, jn
-      !!----------------------------------------------------------------------
-      !
-      !sxM is the first point (in the global domain) needed to compute the north-fold for the current process
-      sxM = jpiglo - nimpp - jpi + 1
-      !dxM is the last point (in the global domain) needed to compute the north-fold for the current process
-      dxM = jpiglo - nimpp + 2
-      !
-      ! loop over the other north-fold processes to find the processes
-      ! managing the points belonging to the sxT-dxT range
-      !
-      nsndto = 0
-      DO jn = 1, jpni
+      !!    3.0  ! 2021-09 complete rewrite using informations from gather north fold
+      !!----------------------------------------------------------------------
+      LOGICAL, INTENT(in   ) ::   ldwrtlay   ! true if additional prints in layout.dat
+      INTEGER, INTENT(in   ) ::   knum       ! layout.dat unit
+      !
+      REAL(wp), DIMENSION(jpi,jpj,2,4) ::   zinfo
+      INTEGER , DIMENSION(10) ::   irknei ! too many elements but safe...
+      INTEGER                 ::   ji, jj, jg, jn   ! dummy loop indices
+      LOGICAL                 ::   lnew
+      !!----------------------------------------------------------------------
+      !
+      IF (lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '   ==>>>   North fold boundary prepared for jpni >1'
+      ENDIF
+      !
+      CALL mpp_ini_northgather   ! we need to init the nfd with gathering in all cases as it is used to define the no-gather case
+      !
+      IF(ldwrtlay) THEN      ! additional prints in layout.dat
+         WRITE(knum,*)
+         WRITE(knum,*)
+         WRITE(knum,*) 'Number of subdomains located along the north fold : ', ndim_rank_north
+         WRITE(knum,*) 'Rank of the subdomains located along the north fold : ', ndim_rank_north
+         DO jn = 1, ndim_rank_north, 5
+            WRITE(knum,*) nrank_north( jn:MINVAL( (/jn+4,ndim_rank_north/) ) )
+         END DO
+      ENDIF
+      
+      nfd_nbnei = 0   ! defaul def (useless?)
+      IF( ln_nnogather ) THEN
          !
-         sxT = nfimpp(jn)                    ! sxT = 1st  point (in the global domain) of the jn process
-         dxT = nfimpp(jn) + nfjpi(jn) - 1    ! dxT = last point (in the global domain) of the jn process
+         ! Use the "gather nfd" to know how to do the nfd: for ji point, which process send data from which of its ji-index?
+         ! Note that nfd is perfectly symetric: I receive data from X <=> I send data to X  (-> no deadlock)
          !
-         IF    ( sxT < sxM  .AND.  sxM < dxT ) THEN
-            nsndto          = nsndto + 1
-            isendto(nsndto) = jn
-         ELSEIF( sxM <= sxT  .AND.  dxM >= dxT ) THEN
-            nsndto          = nsndto + 1
-            isendto(nsndto) = jn
-         ELSEIF( dxM <  dxT  .AND.  sxT <  dxM ) THEN
-            nsndto          = nsndto + 1
-            isendto(nsndto) = jn
-         ENDIF
+         zinfo(:,:,:,:) = HUGE(0._wp)   ! default def to make sure we don't use the halos
+         DO jg = 1, 4   ! grid type: T, U, V, F
+            DO jj = nn_hls+1, jpj-nn_hls                ! inner domain (warning do_loop_substitute not yet defined)
+               DO ji = nn_hls+1, jpi-nn_hls             ! inner domain (warning do_loop_substitute not yet defined)
+                  zinfo(ji,jj,1,jg) = REAL(narea, wp)   ! mpi_rank + 1 (as default lbc_lnk fill is 0
+                  zinfo(ji,jj,2,jg) = REAL(ji, wp)      ! ji of this proc
+               END DO
+            END DO
+         END DO
          !
-      END DO
+         ln_nnogather = .FALSE.   ! force "classical" North pole folding to fill all halos -> should be no more HUGE values...
+         CALL lbc_lnk( 'mppini', zinfo(:,:,:,1), 'T', 1._wp )   ! Do 4 calls instead of 1 to save memory as the nogather version
+         CALL lbc_lnk( 'mppini', zinfo(:,:,:,2), 'U', 1._wp )   ! creates buffer arrays with jpiglo as the first dimension
+         CALL lbc_lnk( 'mppini', zinfo(:,:,:,3), 'V', 1._wp )   ! 
+         CALL lbc_lnk( 'mppini', zinfo(:,:,:,4), 'F', 1._wp )   ! 
+         ln_nnogather = .TRUE.
+         
+         IF( l_IdoNFold ) THEN   ! only the procs involed in the NFD must take care of this
+            
+            ALLOCATE( nfd_rksnd(jpi,4), nfd_jisnd(jpi,4) )          ! neighbour rand and remote ji-index for each grid (T, U, V, F)
+            nfd_rksnd(:,:) = NINT( zinfo(:, jpj, 1, :) ) - 1        ! neighbour MPI rank
+            nfd_jisnd(:,:) = NINT( zinfo(:, jpj, 2, :) ) - nn_hls   ! neighbour ji index (shifted as we don't send the halos)
+            WHERE( nfd_rksnd == -1 )   nfd_jisnd = 1                ! use ji=1 if no neighbour, see mpp_nfd_generic.h90
+            
+            nfd_nbnei = 1                ! Number of neighbour sending data for the nfd. We have at least 1 neighbour!
+            irknei(1) = nfd_rksnd(1,1)   ! which is the 1st one (I can be neighbour of myself, exclude land-proc are also ok)
+            DO jg = 1, 4
+               DO ji = 1, jpi     ! we must be able to fill the full line including halos
+                  lnew = .TRUE.   ! new neighbour?
+                  DO jn = 1, nfd_nbnei
+                     IF( irknei(jn) == nfd_rksnd(ji,jg) )   lnew = .FALSE.   ! already found
+                  END DO
+                  IF( lnew ) THEN
+                     jn = nfd_nbnei + 1
+                     nfd_nbnei = jn
+                     irknei(jn) = nfd_rksnd(ji,jg)
+                  ENDIF
+               END DO
+            END DO
+            
+            ALLOCATE( nfd_rknei(nfd_nbnei) )
+            nfd_rknei(:) = irknei(1:nfd_nbnei)
+            ! re-number nfd_rksnd according to the indexes of nfd_rknei
+            DO jn = 1, nfd_nbnei
+               WHERE( nfd_rksnd == nfd_rknei(jn) )   nfd_rksnd = jn
+            END DO
+            
+            IF( ldwrtlay ) THEN
+               WRITE(knum,*)
+               WRITE(knum,*) 'north fold exchanges with explicit point-to-point messaging :'
+               WRITE(knum,*) '   number of neighbours for the NF: nfd_nbnei : ', nfd_nbnei
+               IF( nfd_nbnei > 0 )   WRITE(knum,*) '   neighbours MPI ranks                       : ', nfd_rknei
+            ENDIF
+            
+         ENDIF   ! l_IdoNFold
+         !
+      ENDIF   ! ln_nnogather
       !
    END SUBROUTINE init_nfdcom
@@ -1342 +1391 @@
    END SUBROUTINE init_doloop
 
+   
+   SUBROUTINE init_locglo
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE init_locglo  ***
+      !!
+      !! ** Purpose :   initialization of global domain <--> local domain indices
+      !!
+      !! ** Method  :
+      !!
+      !! ** Action  : - mig , mjg : local  domain indices ==> global domain, including halos, indices
+      !!              - mig0, mjg0: local  domain indices ==> global domain, excluding halos, indices
+      !!              - mi0 , mi1 : global domain indices ==> local  domain indices
+      !!              - mj0 , mj1   (if global point not in the local domain ==> mi0>mi1 and/or mj0>mj1)
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj   ! dummy loop argument
+      !!----------------------------------------------------------------------
+      !
+      ALLOCATE( mig(jpi), mjg(jpj), mig0(jpi), mjg0(jpj) )
+      ALLOCATE( mi0(jpiglo), mi1(jpiglo), mj0(jpjglo), mj1(jpjglo) )
+      !
+      DO ji = 1, jpi                 ! local domain indices ==> global domain indices, including halos
+        mig(ji) = ji + nimpp - 1
+      END DO
+      DO jj = 1, jpj
+        mjg(jj) = jj + njmpp - 1
+      END DO
+      !                              ! local domain indices ==> global domain indices, excluding halos
+      !
+      mig0(:) = mig(:) - nn_hls
+      mjg0(:) = mjg(:) - nn_hls
+      !                              ! global domain, including halos, indices ==> local domain indices
+      !                                   ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the
+      !                                   ! local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain.
+      DO ji = 1, jpiglo
+        mi0(ji) = MAX( 1 , MIN( ji - nimpp + 1, jpi+1 ) )
+        mi1(ji) = MAX( 0 , MIN( ji - nimpp + 1, jpi   ) )
+      END DO
+      DO jj = 1, jpjglo
+        mj0(jj) = MAX( 1 , MIN( jj - njmpp + 1, jpj+1 ) )
+        mj1(jj) = MAX( 0 , MIN( jj - njmpp + 1, jpj   ) )
+      END DO
+      !
+   END SUBROUTINE init_locglo
+   
    !!======================================================================
 END MODULE mppini

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/SBC/sbcmod.F90

@@ -380 +380 @@
       REAL(wp) ::     zthscl        ! wd  tanh scale
       REAL(wp), DIMENSION(jpi,jpj) ::  zwdht, zwght  ! wd dep over wd limit, wgt
+      REAL(wp), DIMENSION(jpi,jpj) ::  z2d           ! temporary array used for iom_put
 
       !!---------------------------------------------------------------------
@@ -567 +568 @@
       !                                                ! ---------------------------------------- !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-         CALL iom_put( "empmr"  , emp   - rnf )                ! upward water flux
-         CALL iom_put( "empbmr" , emp_b - rnf )                ! before upward water flux ( needed to recalculate the time evolution of ssh in offline )
+         IF( iom_use("empmr") ) THEN
+            DO_2D( 0, 0, 0, 0 )
+               z2d(ji,jj) =  emp(ji,jj) - rnf(ji,jj)
+            END_2D
+            CALL iom_put( "empmr"  , z2d      )                ! upward water flux
+         ENDIF
+         IF( iom_use("empbmr") ) THEN
+            DO_2D( 0, 0, 0, 0 )
+               z2d(ji,jj) =  emp_b(ji,jj) - rnf(ji,jj)
+            END_2D
+            CALL iom_put( "empbmr" , z2d      )                ! before upward water flux ( needed to recalculate the time evolution of ssh in offline )
+         ENDIF
          CALL iom_put( "saltflx", sfx         )                ! downward salt flux (includes virtual salt flux beneath ice in linear free surface case)
          CALL iom_put( "fmmflx" , fmmflx      )                ! Freezing-melting water flux
-         CALL iom_put( "qt"     , qns + qsr   )                ! total heat flux
+         IF( iom_use("qt") ) THEN
+            DO_2D( 0, 0, 0, 0 )
+               z2d(ji,jj) =  qns(ji,jj) + qsr(ji,jj)
+            END_2D
+            CALL iom_put( "qt"  , z2d         )                ! total heat flux
+         ENDIF
          CALL iom_put( "qns"    , qns         )                ! solar heat flux
-         CALL iom_put( "qsr"    ,       qsr   )                ! solar heat flux
+         CALL iom_put( "qsr"    , qsr         )                ! solar heat flux
          IF( nn_ice > 0 .OR. ll_opa )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction
          CALL iom_put( "taum"   , taum        )                ! wind stress module

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/SBC/sbcrnf.F90

@@ -519 +519 @@
       !
       cl_rnfile = TRIM( cn_dir )//TRIM( sn_cnf%clname )
-      IF( .NOT. sn_cnf%ln_clim ) THEN   ;   WRITE(cl_rnfile, '(a,"_y",i4)' ) TRIM( cl_rnfile ), nyear    ! add year
-         IF( sn_cnf%clftyp == 'monthly' )   WRITE(cl_rnfile, '(a,"m",i2)'  ) TRIM( cl_rnfile ), nmonth   ! add month
+      IF( .NOT. sn_cnf%ln_clim ) THEN   ;   WRITE(cl_rnfile, '(a,"_y",i4.4)' ) TRIM( cl_rnfile ), nyear    ! add year
+         IF( sn_cnf%clftyp == 'monthly' )   WRITE(cl_rnfile, '(a,"m" ,i2.2)' ) TRIM( cl_rnfile ), nmonth   ! add month
       ENDIF
       !

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/SBC/sbcssm.F90

@@ -65 +65 @@
       zts(:,:,jp_sal) = ts(:,:,1,jp_sal,Kmm)
       !
-      IF( nn_fsbc == 1 ) THEN                             !   Instantaneous surface fields        !
+         !                                                ! ---------------------------------------- !
+      IF( nn_fsbc == 1 ) THEN                             !      Instantaneous surface fields        !
          !                                                ! ---------------------------------------- !
          ssu_m(:,:) = uu(:,:,1,Kbb)
@@ -92 +93 @@
             ssu_m(:,:) = zcoef * uu(:,:,1,Kbb)
             ssv_m(:,:) = zcoef * vv(:,:,1,Kbb)
-            IF( l_useCT )  THEN    ;   sst_m(:,:) = zcoef * eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
+            IF( l_useCT   )  THEN   ;   sst_m(:,:) = zcoef * eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
             ELSE                    ;   sst_m(:,:) = zcoef * zts(:,:,jp_tem)
             ENDIF

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/SBC/sbcwave.F90

@@ -516 +516 @@
             !
             ! 3. Wave number (only needed for Qiao parametrisation, ln_zdfswm=T)
-            IF( ln_zdfswm .AND. .NOT. cpl_wnum ) THEN
+            IF( .NOT. cpl_wnum ) THEN
                ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wn structure' )

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/TRA/eosbn2.F90

@@ -1081 +1081 @@
          !
          zt  = ctmp   (ji,jj) * z1_T0
-         zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 )
+         zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * z1_S0 )
          ztm = tmask(ji,jj,1)
          !

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/TRA/traadv_mus.F90

@@ -197 +197 @@
          zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions
          zwx(:,:,jpk) = 0._wp
-         DO_3D( 1, 1, 1, 1, 2, jpkm1 )                ! interior values
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! interior values
             zwx(ji,jj,jk) = tmask(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
          END_3D
          !                                !-- Slopes of tracer
          zslpx(:,:,1) = 0._wp                   ! surface values
-         DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             zslpx(ji,jj,jk) =                        ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  &
                &            * (  0.25 + SIGN( 0.25_wp, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  )
          END_3D
-         DO_3D( 1, 1, 1, 1, 2, jpkm1 )    !-- Slopes limitation
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )    !-- Slopes limitation
             zslpx(ji,jj,jk) = SIGN( 1.0_wp, zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
                &                                                     2.*ABS( zwx  (ji,jj,jk+1) ),   &
@@ -221 +221 @@
          IF( ln_linssh ) THEN                   ! top values, linear free surface only
             IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean)
-               DO_2D( 1, 1, 1, 1 )
+               DO_2D( 0, 0, 0, 0 )
                   zwx(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)
                END_2D
             ELSE                                      ! no cavities: only at the ocean surface
-               DO_2D( 1, 1, 1, 1 )
+               DO_2D( 0, 0, 0, 0 )
                   zwx(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb)
                END_2D

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/USR/usrdef_sbc.F90

@@ -110 +110 @@
       ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K)
       zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )   ! emp and rnf used in sshwzv over the whole domain
          ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
          ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
@@ -136 +136 @@
 
       ! freshwater (mass flux) and update of qns with heat content of emp
-      emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average)
-      sfx (:,:) = 0.0_wp                                   ! no salt flux
-      qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )   ! emp used in sshwzv over the whole domain
+         emp (ji,jj) = emp(ji,jj) - zsumemp * tmask(ji,jj,1)          ! freshwater flux (=0 in domain average)
+         sfx (ji,jj) = 0.0_wp                                         ! no salt flux
+         qns (ji,jj) = qns(ji,jj) - emp(ji,jj) * sst_m(ji,jj) * rcp   ! evap and precip are at SST
+      END_2D
 
 
@@ -181 +183 @@
          wndm(ji,jj) = SQRT( zmod * zcoef )
       END_2D
-      CALL lbc_lnk( 'usrdef_sbc', taum(:,:), 'T', 1.0_wp , wndm(:,:), 'T', 1.0_wp )
 
       ! ---------------------------------- !

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/ZDF/zdfevd.F90

@@ -128 +128 @@
          zavt_evd(ji,jj,jk) = p_avt(ji,jj,jk) - zavt_evd(ji,jj,jk)   ! change in avt due to evd
       END_3D
+      !
+      IF( l_trdtra ) CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_evd, zavt_evd )
+      !
       IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
          CALL iom_put( "avt_evd", zavt_evd )              ! output this change
          DEALLOCATE( zavt_evd )
       ENDIF
-      IF( l_trdtra ) CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_evd, zavt_evd )
       !
    END SUBROUTINE zdf_evd

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/ZDF/zdfgls.F90

@@ -1246 +1246 @@
             !
             IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist
-               CALL iom_get( numror, jpdom_auto, 'en'    , en     )
+               CALL iom_get( numror, jpdom_auto, 'en'    , en   , kfill = jpfillcopy  )   ! we devide by en -> must be != 0.
                CALL iom_get( numror, jpdom_auto, 'avt_k' , avt_k  )
                CALL iom_get( numror, jpdom_auto, 'avm_k' , avm_k  )
-               CALL iom_get( numror, jpdom_auto, 'hmxl_n', hmxl_n )
+               CALL iom_get( numror, jpdom_auto, 'hmxl_n', hmxl_n, kfill = jpfillcopy )   ! we devide by hmxl_n -> must be != 0.
             ELSE
                IF(lwp) WRITE(numout,*)

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/ZDF/zdfmxl.F90

@@ -26 +26 @@
    PRIVATE
 
-   PUBLIC   zdf_mxl, zdf_mxl_turb   ! called by zdfphy.F90
+   PUBLIC   zdf_mxl, zdf_mxl_turb, zdf_mxl_alloc   ! called by zdfphy.F90
 
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by LDF, ZDF, TRD, TOP)
@@ -86 +86 @@
             IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
             IF(lwp) WRITE(numout,*) '~~~~~~~ '
-            !                             ! allocate zdfmxl arrays
-            IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
          ENDIF
       ENDIF
@@ -111 +109 @@
       !
       IF( .NOT.l_offline .AND. iom_use("mldr10_1") ) THEN
-         IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness
-         ELSE                  ;  CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth
-         END IF
+         IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN         ! Do only on the last tile
+            IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness
+            ELSE                  ;  CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth
+            END IF
+         ENDIF
       ENDIF
       !
@@ -144 +144 @@
       ! w-level of the turbocline and mixing layer (iom_use)
       imld(:,:) = mbkt(A2D(nn_hls)) + 1                ! Initialization to the number of w ocean point
-      DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, nlb10, -1 )   ! from the bottom to nlb10
+      DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )   ! from the bottom to nlb10
          IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline
       END_3D
       ! depth of the mixing layer
-      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D_OVR( 1, 1, 1, 1 )
          iik = imld(ji,jj)
          hmld (ji,jj) = gdepw(ji,jj,iik  ,Kmm) * ssmask(ji,jj)    ! Turbocline depth
@@ -154 +154 @@
       !
       IF( .NOT.l_offline .AND. iom_use("mldkz5") ) THEN
-         IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness
-         ELSE                  ;  CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth
-         END IF
+         IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN         ! Do only on the last tile
+            IF( ln_isfcav ) THEN  ;  CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness
+            ELSE                  ;  CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth
+            END IF
+         ENDIF
       ENDIF
       !

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/ZDF/zdfphy.F90

@@ -146 +146 @@
          wi    (:,:,:) = 0._wp
       ENDIF
+      !                                  ! Initialise zdf_mxl arrays (only hmld as not set everywhere when nn_hls > 1)
+      IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl : unable to allocate arrays' )
+      hmld(:,:) = 0._wp
       !                          !==  Background eddy viscosity and diffusivity  ==!
       IF( nn_avb == 0 ) THEN             ! Define avmb, avtb from namelist parameter
@@ -378 +381 @@
       IF( iom_use('avt_k') .OR. iom_use('avm_k') .OR. iom_use('eshear_k') .OR. iom_use('estrat_k') ) THEN
          IF( l_zdfsh2 ) THEN
-            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
-               zsh2(ji,jj,1  ) = 0._wp
-               zsh2(ji,jj,jpk) = 0._wp
-            END_2D
             CALL iom_put( 'avt_k'   ,   avt_k       * wmask )
             CALL iom_put( 'avm_k'   ,   avm_k       * wmask )

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/ZDF/zdfric.F90

@@ -157 +157 @@
       !
       !                       !==  avm and avt = F(Richardson number)  ==!
-      DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, jpkm1 )       ! coefficient = F(richardson number) (avm-weighted Ri)
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! coefficient = F(richardson number) (avm-weighted Ri)
          zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) )  )
          zav   = rn_avmri * zcfRi**nn_ric

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/ZDF/zdfsh2.F90

@@ -96 +96 @@
          END_2D
       END DO
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! set p_sh2 to 0 at the surface and bottom for output purpose
+         p_sh2(ji,jj,1)   = 0._wp
+         p_sh2(ji,jj,jpk) = 0._wp
+      END_2D
       !
    END SUBROUTINE zdf_sh2

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/lib_fortran.F90

@@ -26 +26 @@
    PRIVATE
 
-   PUBLIC   glob_sum      ! used in many places (masked with tmask_i)
-   PUBLIC   glob_sum_full ! used in many places (masked with tmask_h, ie only over the halos)
+   PUBLIC   glob_sum      ! used in many places (masked with tmask_i = ssmask * tmask_h)
+   PUBLIC   glob_sum_full ! used in many places (masked with tmask_h, excluding all duplicated points halos+periodicity)
    PUBLIC   local_sum     ! used in trcrad, local operation before glob_sum_delay
    PUBLIC   sum3x3        ! used in trcrad, do a sum over 3x3 boxes
@@ -331 +331 @@
       INTEGER     ::   ji , jj , jk     ! dummy loop indices
       INTEGER     ::   ipi, ipj, ipk    ! dimensions
+      INTEGER     ::   iis, iie, ijs, ije   ! loop start and end
       !!-----------------------------------------------------------------------
       !
@@ -337 +338 @@
       ipk = SIZE(ptab,3)   ! 3rd dimension
       !
+      IF( ipi == jpi .AND. ipj == jpj ) THEN   ! do 2D loop only over the inner domain (-> avoid to use undefined values)
+         iis = Nis0   ;   iie = Nie0
+         ijs = Njs0   ;   ije = Nje0
+      ELSE                                     ! I think we are never in this case...
+         iis = 1   ;   iie = jpi
+         ijs = 1   ;   ije = jpj
+      ENDIF
+      !
       ALLOCATE( ctmp(ipk) )
       !
       DO jk = 1, ipk
          ctmp(jk) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
-         DO jj = 1, ipj
-            DO ji = 1, ipi
+         DO jj = ijs, ije
+            DO ji = iis, iie
                ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
                CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
@@ -366 +375 @@
       INTEGER     ::   ji , jj , jk , jl     ! dummy loop indices
       INTEGER     ::   ipi, ipj, ipk, ipl    ! dimensions
+      INTEGER     ::   iis, iie, ijs, ije    ! loop start and end
       !!-----------------------------------------------------------------------
       !
@@ -373 +383 @@
       ipl = SIZE(ptab,4)   ! 4th dimension
       !
+      IF( ipi == jpi .AND. ipj == jpj ) THEN   ! do 2D loop only over the inner domain (-> avoid to use undefined values)
+         iis = Nis0   ;   iie = Nie0
+         ijs = Njs0   ;   ije = Nje0
+      ELSE                                     ! I think we are never in this case...
+         iis = 1   ;   iie = jpi
+         ijs = 1   ;   ije = jpj
+      ENDIF
+      !
       ALLOCATE( ctmp(ipl) )
       !
@@ -378 +396 @@
          ctmp(jl) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
          DO jk = 1, ipk
-            DO jj = 1, ipj
-               DO ji = 1, ipi
+            DO jj = ijs, ije
+               DO ji = iis, iie
                   ztmp =  ptab(ji,jj,jk,jl) * tmask_i(ji,jj)
                   CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )
@@ -404 +422 @@
       INTEGER     ::   ji , jj , jk     ! dummy loop indices
       INTEGER     ::   ipi, ipj, ipk    ! dimensions
+      INTEGER     ::   iis, iie, ijs, ije   ! loop start and end
       !!-----------------------------------------------------------------------
       !
@@ -410 +429 @@
       ipk = SIZE(ptab,3)   ! 3rd dimension
       !
+      IF( ipi == jpi .AND. ipj == jpj ) THEN   ! do 2D loop only over the inner domain (-> avoid to use undefined values)
+         iis = Nis0   ;   iie = Nie0
+         ijs = Njs0   ;   ije = Nje0
+      ELSE                                     ! I think we are never in this case...
+         iis = 1   ;   iie = jpi
+         ijs = 1   ;   ije = jpj
+      ENDIF
+      !
       ALLOCATE( ctmp(ipk) )
       !
       DO jk = 1, ipk
          ctmp(jk) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
-         DO jj = 1, ipj
-            DO ji = 1, ipi
+         DO jj = ijs, ije
+            DO ji = iis, iie
                ztmp =  ptab(ji,jj,jk) * tmask_h(ji,jj)
                CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
@@ -439 +466 @@
       INTEGER     ::   ji , jj , jk , jl     ! dummy loop indices
       INTEGER     ::   ipi, ipj, ipk, ipl    ! dimensions
+      INTEGER     ::   iis, iie, ijs, ije    ! loop start and end
       !!-----------------------------------------------------------------------
       !
@@ -446 +474 @@
       ipl = SIZE(ptab,4)   ! 4th dimension
       !
+      IF( ipi == jpi .AND. ipj == jpj ) THEN   ! do 2D loop only over the inner domain (-> avoid to use undefined values)
+         iis = Nis0   ;   iie = Nie0
+         ijs = Njs0   ;   ije = Nje0
+      ELSE                                     ! I think we are never in this case...
+         iis = 1   ;   iie = jpi
+         ijs = 1   ;   ije = jpj
+      ENDIF
+      !
       ALLOCATE( ctmp(ipl) )
       !
@@ -451 +487 @@
          ctmp(jl) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
          DO jk = 1, ipk
-            DO jj = 1, ipj
-               DO ji = 1, ipi
+            DO jj = ijs, ije
+               DO ji = iis, iie
                   ztmp =  ptab(ji,jj,jk,jl) * tmask_h(ji,jj)
                   CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/lib_fortran_generic.h90

@@ -38 +38 @@
       !!-----------------------------------------------------------------------