New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 13159 for NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP+ – NEMO

Ignore:
Timestamp:
2020-06-26T10:26:32+02:00 (4 years ago)
Author:
gsamson
Message:

merge trunk@r13136 into ASINTER-06 branch; pass all SETTE tests; results identical to trunk@r13136; ticket #2419

Location:
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement

    • Property svn:externals
      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@HEAD         sette 
         10^/utils/CI/sette@12931        sette 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP+/EXPREF/file_def_nemo-oce.xml

    r11889 r13159  
    2121      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE.">  <!-- 5d files -->   
    2222 
    23       <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files --> 
    24    <file id="file1" output_freq="1mo" name_suffix="_grid_T" description="ocean T grid variables" > 
    25      <field field_ref="toce"         name="votemper"  /> 
    26      <field field_ref="soce"         name="vosaline"  /> 
    27      <field field_ref="ssh"          name="sossheig"  /> 
     23   <file id="file1" output_freq="5d" name_suffix="_grid_T" description="ocean T grid variables" > 
     24     <field field_ref="toce"         name="votemper"  operation="average" freq_op="5d" > @toce_e3t / @e3t </field> 
     25     <field field_ref="soce"         name="vosaline"  operation="average" freq_op="5d" > @soce_e3t / @e3t </field> 
     26     <field field_ref="ssh"          name="sossheig" /> 
    2827          <!-- variable for ice shelf --> 
    29           <field field_ref="fwfisf_cav"       name="sowflisf"  /> 
    30           <field field_ref="isfgammat"    name="sogammat"  /> 
    31           <field field_ref="isfgammas"    name="sogammas"  /> 
     28          <field field_ref="fwfisf_cav"  name="sowflisf"  /> 
     29          <field field_ref="isfgammat"   name="sogammat"  /> 
     30          <field field_ref="isfgammas"   name="sogammas"  /> 
    3231          <field field_ref="ttbl_cav"    name="ttbl"  /> 
    33           <field field_ref="stbl"    name="stbl"  /> 
    34           <field field_ref="utbl"    name="utbl"  /> 
    35           <field field_ref="vtbl"    name="vtbl"  /> 
     32          <field field_ref="stbl"        name="stbl"  /> 
     33          <field field_ref="utbl"        name="utbl"  /> 
     34          <field field_ref="vtbl"        name="vtbl"  /> 
    3635        </file> 
    37    <file id="file2" output_freq="1mo" name_suffix="_grid_U" description="ocean U grid variables" > 
    38           <field field_ref="uoce"         name="vozocrtx" /> 
     36   <file id="file2" output_freq="5d" name_suffix="_grid_U" description="ocean U grid variables" > 
     37          <field field_ref="uoce"         name="vozocrtx" operation="average" freq_op="5d" > @uoce_e3u / @e3u </field> /> 
    3938        </file> 
    40    <file id="file3" output_freq="1mo" name_suffix="_grid_V" description="ocean V grid variables" > 
    41           <field field_ref="voce"         name="vomecrty" />  
     39   <file id="file3" output_freq="5d" name_suffix="_grid_V" description="ocean V grid variables" > 
     40          <field field_ref="voce"         name="vomecrty" operation="average" freq_op="5d" > @voce_e3v / @e3v </field> />  
    4241        </file> 
    4342      </file_group> 
     43 
     44      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files --> 
    4445      <file_group id="2m" output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files --> 
    4546      <file_group id="3m" output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files --> 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP+/EXPREF/namelist_cfg

    r12489 r13159  
    114114 
    115115   ln_usr      = .true.   !  user defined formulation                  (T => check usrdef_sbc) 
    116    nn_fwb      = 1 
     116   nn_fwb      = 4 
    117117/ 
    118118!----------------------------------------------------------------------- 
     
    308308&nameos        !   ocean Equation Of Seawater                           (default: NO selection) 
    309309!----------------------------------------------------------------------- 
    310    ln_teos10   = .false.         !  = Use TEOS-10 
    311    ln_eos80    = .false.         !  = Use EOS80 
    312    ln_leos     = .true.          !  = Use S-EOS (simplified Eq.) 
     310   ln_leos     = .true.          !  = Use L-EOS (linear Eq.) 
    313311                                 ! 
    314312   !                     ! S-EOS coefficients (ln_seos=T): 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP+/MY_SRC/dtatsd.F90

    r12077 r13159  
    3636   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsddmp ! structure of input SST (file informations, fields read) 
    3737 
     38   !! * Substitutions 
     39#  include "do_loop_substitute.h90" 
    3840   !!---------------------------------------------------------------------- 
    3941   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    6769      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0 
    6870      ! 
    69       REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :  
    7071      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901) 
    7172901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtsd in reference namelist' ) 
    72       REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run 
    7373      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 ) 
    7474902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtsd in configuration namelist' ) 
     
    191191         ENDIF 
    192192         ! 
    193          DO jj = 1, jpj                         ! vertical interpolation of T & S 
    194             DO ji = 1, jpi 
    195                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
    196                   zl = gdept_0(ji,jj,jk) 
    197                   IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data 
    198                      ztp(jk) =  ptsd(ji,jj,1    ,jp_tem) 
    199                      zsp(jk) =  ptsd(ji,jj,1    ,jp_sal) 
    200                   ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data 
    201                      ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem) 
    202                      zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal) 
    203                   ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
    204                      DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
    205                         IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
    206                            zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
    207                            ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi  
    208                            zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi 
    209                         ENDIF 
    210                      END DO 
    211                   ENDIF 
    212                END DO 
    213                DO jk = 1, jpkm1 
    214                   ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
    215                   ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk) 
    216                END DO 
    217                ptsd(ji,jj,jpk,jp_tem) = 0._wp 
    218                ptsd(ji,jj,jpk,jp_sal) = 0._wp 
     193         DO_2D_11_11 
     194            DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
     195               zl = gdept_0(ji,jj,jk) 
     196               IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data 
     197                  ztp(jk) =  ptsd(ji,jj,1    ,jp_tem) 
     198                  zsp(jk) =  ptsd(ji,jj,1    ,jp_sal) 
     199               ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data 
     200                  ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem) 
     201                  zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal) 
     202               ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
     203                  DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
     204                     IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
     205                        zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
     206                        ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi  
     207                        zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi 
     208                     ENDIF 
     209                  END DO 
     210               ENDIF 
    219211            END DO 
    220          END DO 
     212            DO jk = 1, jpkm1 
     213               ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
     214               ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk) 
     215            END DO 
     216            ptsd(ji,jj,jpk,jp_tem) = 0._wp 
     217            ptsd(ji,jj,jpk,jp_sal) = 0._wp 
     218         END_2D 
    221219         !  
    222220      ELSE                                !==   z- or zps- coordinate   ==! 
     
    226224         ! 
    227225         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level 
    228             DO jj = 1, jpj 
    229                DO ji = 1, jpi 
    230                   ik = mbkt(ji,jj)  
    231                   IF( ik > 1 ) THEN 
    232                      zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
    233                      ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem) 
    234                      ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal) 
    235                   ENDIF 
    236                   ik = mikt(ji,jj) 
    237                   IF( ik > 1 ) THEN 
    238                      zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )  
    239                      ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem) 
    240                      ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal) 
    241                   END IF 
    242                END DO 
    243             END DO 
     226            DO_2D_11_11 
     227               ik = mbkt(ji,jj)  
     228               IF( ik > 1 ) THEN 
     229                  zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
     230                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem) 
     231                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal) 
     232               ENDIF 
     233               ik = mikt(ji,jj) 
     234               IF( ik > 1 ) THEN 
     235                  zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )  
     236                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem) 
     237                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal) 
     238               END IF 
     239            END_2D 
    244240         ENDIF 
    245241         ! 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP+/MY_SRC/eosbn2.F90

    r12489 r13159  
    180180   REAL(wp) ::   BPE002 
    181181 
     182   !! * Substitutions 
     183#  include "do_loop_substitute.h90" 
    182184   !!---------------------------------------------------------------------- 
    183185   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    241243      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    242244         ! 
    243          DO jk = 1, jpkm1 
    244             DO jj = 1, jpj 
    245                DO ji = 1, jpi 
    246                   ! 
    247                   zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    248                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    249                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    250                   ztm = tmask(ji,jj,jk)                                         ! tmask 
     245         DO_3D_11_11( 1, jpkm1 ) 
     246            ! 
     247            zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     248            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     249            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     250            ztm = tmask(ji,jj,jk)                                         ! tmask 
     251            ! 
     252            zn3 = EOS013*zt   & 
     253               &   + EOS103*zs+EOS003 
     254               ! 
     255            zn2 = (EOS022*zt   & 
     256               &   + EOS112*zs+EOS012)*zt   & 
     257               &   + (EOS202*zs+EOS102)*zs+EOS002 
     258               ! 
     259            zn1 = (((EOS041*zt   & 
     260               &   + EOS131*zs+EOS031)*zt   & 
     261               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     262               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     263               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     264               ! 
     265            zn0 = (((((EOS060*zt   & 
     266               &   + EOS150*zs+EOS050)*zt   & 
     267               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     268               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     269               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     270               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     271               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     272               ! 
     273            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     274            ! 
     275            prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked) 
     276            ! 
     277         END_3D 
     278         ! 
     279      CASE( np_seos )                !==  simplified EOS  ==! 
     280         ! 
     281         DO_3D_11_11( 1, jpkm1 ) 
     282            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
     283            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
     284            zh  = pdep (ji,jj,jk) 
     285            ztm = tmask(ji,jj,jk) 
     286            ! 
     287            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
     288               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
     289               &  - rn_nu * zt * zs 
     290               !                                  
     291            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked) 
     292         END_3D 
     293         ! 
     294      CASE( np_leos )                !==  linear ISOMIP EOS  ==! 
     295         ! 
     296         DO_3D_11_11( 1, jpkm1 ) 
     297            zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp) 
     298            zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp 
     299            zh  = pdep (ji,jj,jk) 
     300            ztm = tmask(ji,jj,jk) 
     301            ! 
     302            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
     303            !                                  
     304            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked) 
     305         END_3D 
     306         ! 
     307      END SELECT 
     308      ! 
     309      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk ) 
     310      ! 
     311      IF( ln_timing )   CALL timing_stop('eos-insitu') 
     312      ! 
     313   END SUBROUTINE eos_insitu 
     314 
     315 
     316   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 
     317      !!---------------------------------------------------------------------- 
     318      !!                  ***  ROUTINE eos_insitu_pot  *** 
     319      !! 
     320      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the 
     321      !!      potential volumic mass (Kg/m3) from potential temperature and 
     322      !!      salinity fields using an equation of state selected in the 
     323      !!     namelist. 
     324      !! 
     325      !! ** Action  : - prd  , the in situ density (no units) 
     326      !!              - prhop, the potential volumic mass (Kg/m3) 
     327      !! 
     328      !!---------------------------------------------------------------------- 
     329      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
     330      !                                                                ! 2 : salinity               [psu] 
     331      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-] 
     332      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
     333      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
     334      ! 
     335      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices 
     336      INTEGER  ::   jdof 
     337      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars 
     338      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      - 
     339      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors 
     340      !!---------------------------------------------------------------------- 
     341      ! 
     342      IF( ln_timing )   CALL timing_start('eos-pot') 
     343      ! 
     344      SELECT CASE ( neos ) 
     345      ! 
     346      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     347         ! 
     348         ! Stochastic equation of state 
     349         IF ( ln_sto_eos ) THEN 
     350            ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
     351            ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
     352            ALLOCATE(zsign(1:2*nn_sto_eos)) 
     353            DO jsmp = 1, 2*nn_sto_eos, 2 
     354              zsign(jsmp)   = 1._wp 
     355              zsign(jsmp+1) = -1._wp 
     356            END DO 
     357            ! 
     358            DO_3D_11_11( 1, jpkm1 ) 
     359               ! 
     360               ! compute density (2*nn_sto_eos) times: 
     361               ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
     362               ! (2) for t-dt, s-ds (with the opposite fluctuation) 
     363               DO jsmp = 1, nn_sto_eos*2 
     364                  jdof   = (jsmp + 1) / 2 
     365                  zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     366                  zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature 
     367                  zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
     368                  zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
     369                  ztm    = tmask(ji,jj,jk)                                         ! tmask 
    251370                  ! 
    252371                  zn3 = EOS013*zt   & 
     
    263382                     &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    264383                     ! 
    265                   zn0 = (((((EOS060*zt   & 
     384                  zn0_sto(jsmp) = (((((EOS060*zt   & 
    266385                     &   + EOS150*zs+EOS050)*zt   & 
    267386                     &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     
    271390                     &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    272391                     ! 
    273                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     392                  zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 
     393               END DO 
     394               ! 
     395               ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
     396               prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
     397               DO jsmp = 1, nn_sto_eos*2 
     398                  prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
    274399                  ! 
    275                   prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked) 
    276                   ! 
     400                  prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked) 
    277401               END DO 
    278             END DO 
    279          END DO 
    280          ! 
    281       CASE( np_seos )                !==  simplified EOS  ==! 
    282          ! 
    283          DO jk = 1, jpkm1 
    284             DO jj = 1, jpj 
    285                DO ji = 1, jpi 
    286                   zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
    287                   zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
    288                   zh  = pdep (ji,jj,jk) 
    289                   ztm = tmask(ji,jj,jk) 
    290                   ! 
    291                   zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
    292                      &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
    293                      &  - rn_nu * zt * zs 
    294                      !                                  
    295                   prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked) 
    296                END DO 
    297             END DO 
    298          END DO 
    299          ! 
    300       CASE( np_leos )                !==  linear ISOMIP EOS  ==! 
    301          ! 
    302          DO jk = 1, jpkm1 
    303             DO jj = 1, jpj 
    304                DO ji = 1, jpi 
    305                   zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp) 
    306                   zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp 
    307                   zh  = pdep (ji,jj,jk) 
    308                   ztm = tmask(ji,jj,jk) 
    309                   ! 
    310                   zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
    311                   !                                  
    312                   prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked) 
    313                END DO 
    314             END DO 
    315          END DO 
    316          ! 
    317       END SELECT 
    318       ! 
    319       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk ) 
    320       ! 
    321       IF( ln_timing )   CALL timing_stop('eos-insitu') 
    322       ! 
    323    END SUBROUTINE eos_insitu 
    324  
    325  
    326    SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 
    327       !!---------------------------------------------------------------------- 
    328       !!                  ***  ROUTINE eos_insitu_pot  *** 
    329       !! 
    330       !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the 
    331       !!      potential volumic mass (Kg/m3) from potential temperature and 
    332       !!      salinity fields using an equation of state selected in the 
    333       !!     namelist. 
    334       !! 
    335       !! ** Action  : - prd  , the in situ density (no units) 
    336       !!              - prhop, the potential volumic mass (Kg/m3) 
    337       !! 
    338       !!---------------------------------------------------------------------- 
    339       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
    340       !                                                                ! 2 : salinity               [psu] 
    341       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-] 
    342       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
    343       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
    344       ! 
    345       INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices 
    346       INTEGER  ::   jdof 
    347       REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars 
    348       REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      - 
    349       REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors 
    350       !!---------------------------------------------------------------------- 
    351       ! 
    352       IF( ln_timing )   CALL timing_start('eos-pot') 
    353       ! 
    354       SELECT CASE ( neos ) 
    355       ! 
    356       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    357          ! 
    358          ! Stochastic equation of state 
    359          IF ( ln_sto_eos ) THEN 
    360             ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
    361             ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
    362             ALLOCATE(zsign(1:2*nn_sto_eos)) 
    363             DO jsmp = 1, 2*nn_sto_eos, 2 
    364               zsign(jsmp)   = 1._wp 
    365               zsign(jsmp+1) = -1._wp 
    366             END DO 
    367             ! 
    368             DO jk = 1, jpkm1 
    369                DO jj = 1, jpj 
    370                   DO ji = 1, jpi 
    371                      ! 
    372                      ! compute density (2*nn_sto_eos) times: 
    373                      ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
    374                      ! (2) for t-dt, s-ds (with the opposite fluctuation) 
    375                      DO jsmp = 1, nn_sto_eos*2 
    376                         jdof   = (jsmp + 1) / 2 
    377                         zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    378                         zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature 
    379                         zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
    380                         zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
    381                         ztm    = tmask(ji,jj,jk)                                         ! tmask 
    382                         ! 
    383                         zn3 = EOS013*zt   & 
    384                            &   + EOS103*zs+EOS003 
    385                            ! 
    386                         zn2 = (EOS022*zt   & 
    387                            &   + EOS112*zs+EOS012)*zt   & 
    388                            &   + (EOS202*zs+EOS102)*zs+EOS002 
    389                            ! 
    390                         zn1 = (((EOS041*zt   & 
    391                            &   + EOS131*zs+EOS031)*zt   & 
    392                            &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
    393                            &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
    394                            &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    395                            ! 
    396                         zn0_sto(jsmp) = (((((EOS060*zt   & 
    397                            &   + EOS150*zs+EOS050)*zt   & 
    398                            &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    399                            &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
    400                            &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
    401                            &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
    402                            &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    403                            ! 
    404                         zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 
    405                      END DO 
    406                      ! 
    407                      ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
    408                      prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
    409                      DO jsmp = 1, nn_sto_eos*2 
    410                         prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
    411                         ! 
    412                         prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked) 
    413                      END DO 
    414                      prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
    415                      prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
    416                   END DO 
    417                END DO 
    418             END DO 
     402               prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
     403               prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
     404            END_3D 
    419405            DEALLOCATE(zn0_sto,zn_sto,zsign) 
    420406         ! Non-stochastic equation of state 
    421407         ELSE 
    422             DO jk = 1, jpkm1 
    423                DO jj = 1, jpj 
    424                   DO ji = 1, jpi 
    425                      ! 
    426                      zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    427                      zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    428                      zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    429                      ztm = tmask(ji,jj,jk)                                         ! tmask 
    430                      ! 
    431                      zn3 = EOS013*zt   & 
    432                         &   + EOS103*zs+EOS003 
    433                         ! 
    434                      zn2 = (EOS022*zt   & 
    435                         &   + EOS112*zs+EOS012)*zt   & 
    436                         &   + (EOS202*zs+EOS102)*zs+EOS002 
    437                         ! 
    438                      zn1 = (((EOS041*zt   & 
    439                         &   + EOS131*zs+EOS031)*zt   & 
    440                         &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
    441                         &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
    442                         &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    443                         ! 
    444                      zn0 = (((((EOS060*zt   & 
    445                         &   + EOS150*zs+EOS050)*zt   & 
    446                         &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    447                         &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
    448                         &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
    449                         &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
    450                         &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    451                         ! 
    452                      zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    453                      ! 
    454                      prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
    455                      ! 
    456                      prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked) 
    457                   END DO 
    458                END DO 
    459             END DO 
    460          ENDIF 
    461           
    462       CASE( np_seos )                !==  simplified EOS  ==! 
    463          ! 
    464          DO jk = 1, jpkm1 
    465             DO jj = 1, jpj 
    466                DO ji = 1, jpi 
    467                   zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
    468                   zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
    469                   zh  = pdep (ji,jj,jk) 
    470                   ztm = tmask(ji,jj,jk) 
    471                   !                                                     ! potential density referenced at the surface 
    472                   zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   & 
    473                      &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   & 
    474                      &  - rn_nu * zt * zs 
    475                   prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 
    476                   !                                                     ! density anomaly (masked) 
    477                   zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 
    478                   prd(ji,jj,jk) = zn * r1_rho0 * ztm 
    479                   ! 
    480                END DO 
    481             END DO 
    482          END DO 
    483          ! 
    484       CASE( np_leos )                !==  linear ISOMIP EOS  ==! 
    485          ! 
    486          DO jk = 1, jpkm1 
    487             DO jj = 1, jpj 
    488                DO ji = 1, jpi 
    489                   zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp) 
    490                   zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp 
    491                   zh  = pdep (ji,jj,jk) 
    492                   ztm = tmask(ji,jj,jk) 
    493                   !                                                     ! potential density referenced at the surface 
    494                   zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
    495                   prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 
    496                   !                                                     ! density anomaly (masked) 
    497                   prd(ji,jj,jk) = zn * r1_rho0 * ztm 
    498                   ! 
    499                END DO 
    500             END DO 
    501          END DO 
    502          ! 
    503       END SELECT 
    504       ! 
    505       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk ) 
    506       ! 
    507       IF( ln_timing )   CALL timing_stop('eos-pot') 
    508       ! 
    509    END SUBROUTINE eos_insitu_pot 
    510  
    511  
    512    SUBROUTINE eos_insitu_2d( pts, pdep, prd ) 
    513       !!---------------------------------------------------------------------- 
    514       !!                  ***  ROUTINE eos_insitu_2d  *** 
    515       !! 
    516       !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from 
    517       !!      potential temperature and salinity using an equation of state 
    518       !!      selected in the nameos namelist. * 2D field case 
    519       !! 
    520       !! ** Action  : - prd , the in situ density (no units) (unmasked) 
    521       !! 
    522       !!---------------------------------------------------------------------- 
    523       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius] 
    524       !                                                           ! 2 : salinity               [psu] 
    525       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m] 
    526       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
    527       ! 
    528       INTEGER  ::   ji, jj, jk                ! dummy loop indices 
    529       REAL(wp) ::   zt , zh , zs              ! local scalars 
    530       REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
    531       !!---------------------------------------------------------------------- 
    532       ! 
    533       IF( ln_timing )   CALL timing_start('eos2d') 
    534       ! 
    535       prd(:,:) = 0._wp 
    536       ! 
    537       SELECT CASE( neos ) 
    538       ! 
    539       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    540          ! 
    541          DO jj = 1, jpjm1 
    542             DO ji = 1, fs_jpim1   ! vector opt. 
    543                ! 
    544                zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
    545                zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
    546                zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     408            DO_3D_11_11( 1, jpkm1 ) 
     409               ! 
     410               zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     411               zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     412               zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     413               ztm = tmask(ji,jj,jk)                                         ! tmask 
    547414               ! 
    548415               zn3 = EOS013*zt   & 
     
    569436               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    570437               ! 
    571                prd(ji,jj) = zn * r1_rho0 - 1._wp               ! unmasked in situ density anomaly 
    572                ! 
    573             END DO 
    574          END DO 
    575          ! 
    576          CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions 
    577          ! 
     438               prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
     439               ! 
     440               prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked) 
     441            END_3D 
     442         ENDIF 
     443          
    578444      CASE( np_seos )                !==  simplified EOS  ==! 
    579445         ! 
    580          DO jj = 1, jpjm1 
    581             DO ji = 1, fs_jpim1   ! vector opt. 
    582                ! 
    583                zt    = pts  (ji,jj,jp_tem)  - 10._wp 
    584                zs    = pts  (ji,jj,jp_sal)  - 35._wp 
    585                zh    = pdep (ji,jj)                         ! depth at the partial step level 
    586                ! 
    587                zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
    588                   &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
    589                   &  - rn_nu * zt * zs 
    590                   ! 
    591                prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly 
    592                ! 
    593             END DO 
    594          END DO 
    595          ! 
    596          CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions 
     446         DO_3D_11_11( 1, jpkm1 ) 
     447            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
     448            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
     449            zh  = pdep (ji,jj,jk) 
     450            ztm = tmask(ji,jj,jk) 
     451            !                                                     ! potential density referenced at the surface 
     452            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   & 
     453               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   & 
     454               &  - rn_nu * zt * zs 
     455            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 
     456            !                                                     ! density anomaly (masked) 
     457            zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 
     458            prd(ji,jj,jk) = zn * r1_rho0 * ztm 
     459            ! 
     460         END_3D 
     461         ! 
     462      CASE( np_leos )                !==  linear ISOMIP EOS  ==! 
     463         ! 
     464         DO_3D_11_11( 1, jpkm1 ) 
     465            zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp) 
     466            zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp 
     467            zh  = pdep (ji,jj,jk) 
     468            ztm = tmask(ji,jj,jk) 
     469            !                                                     ! potential density referenced at the surface 
     470            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
     471            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 
     472            !                                                     ! density anomaly (masked) 
     473            prd(ji,jj,jk) = zn * r1_rho0 * ztm 
     474            ! 
     475         END_3D 
     476         ! 
     477      END SELECT 
     478      ! 
     479      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk ) 
     480      ! 
     481      IF( ln_timing )   CALL timing_stop('eos-pot') 
     482      ! 
     483   END SUBROUTINE eos_insitu_pot 
     484 
     485 
     486   SUBROUTINE eos_insitu_2d( pts, pdep, prd ) 
     487      !!---------------------------------------------------------------------- 
     488      !!                  ***  ROUTINE eos_insitu_2d  *** 
     489      !! 
     490      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from 
     491      !!      potential temperature and salinity using an equation of state 
     492      !!      selected in the nameos namelist. * 2D field case 
     493      !! 
     494      !! ** Action  : - prd , the in situ density (no units) (unmasked) 
     495      !! 
     496      !!---------------------------------------------------------------------- 
     497      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius] 
     498      !                                                           ! 2 : salinity               [psu] 
     499      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m] 
     500      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
     501      ! 
     502      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     503      REAL(wp) ::   zt , zh , zs              ! local scalars 
     504      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     505      !!---------------------------------------------------------------------- 
     506      ! 
     507      IF( ln_timing )   CALL timing_start('eos2d') 
     508      ! 
     509      prd(:,:) = 0._wp 
     510      ! 
     511      SELECT CASE( neos ) 
     512      ! 
     513      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     514         ! 
     515         DO_2D_11_11 
     516            ! 
     517            zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
     518            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
     519            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     520            ! 
     521            zn3 = EOS013*zt   & 
     522               &   + EOS103*zs+EOS003 
     523               ! 
     524            zn2 = (EOS022*zt   & 
     525               &   + EOS112*zs+EOS012)*zt   & 
     526               &   + (EOS202*zs+EOS102)*zs+EOS002 
     527               ! 
     528            zn1 = (((EOS041*zt   & 
     529               &   + EOS131*zs+EOS031)*zt   & 
     530               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     531               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     532               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     533               ! 
     534            zn0 = (((((EOS060*zt   & 
     535               &   + EOS150*zs+EOS050)*zt   & 
     536               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     537               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     538               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     539               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     540               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     541               ! 
     542            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     543            ! 
     544            prd(ji,jj) = zn * r1_rho0 - 1._wp               ! unmasked in situ density anomaly 
     545            ! 
     546         END_2D 
     547         ! 
     548      CASE( np_seos )                !==  simplified EOS  ==! 
     549         ! 
     550         DO_2D_11_11 
     551            ! 
     552            zt    = pts  (ji,jj,jp_tem)  - 10._wp 
     553            zs    = pts  (ji,jj,jp_sal)  - 35._wp 
     554            zh    = pdep (ji,jj)                         ! depth at the partial step level 
     555            ! 
     556            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
     557               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
     558               &  - rn_nu * zt * zs 
     559               ! 
     560            prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly 
     561            ! 
     562         END_2D 
    597563         ! 
    598564      CASE( np_leos )                !==  ISOMIP EOS  ==! 
    599565         ! 
    600          DO jj = 1, jpjm1 
    601             DO ji = 1, fs_jpim1   ! vector opt. 
    602                ! 
    603                zt    = pts  (ji,jj,jp_tem)  - (-1._wp) 
    604                zs    = pts  (ji,jj,jp_sal)  - 34.2_wp 
    605                zh    = pdep (ji,jj)                         ! depth at the partial step level 
    606                ! 
    607                zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
    608                   ! 
    609                prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly 
    610                ! 
    611             END DO 
    612          END DO 
    613          ! 
    614          CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions 
     566         DO_2D_11_11 
     567            ! 
     568            zt    = pts  (ji,jj,jp_tem)  - (-1._wp) 
     569            zs    = pts  (ji,jj,jp_sal)  - 34.2_wp 
     570            zh    = pdep (ji,jj)                         ! depth at the partial step level 
     571            ! 
     572            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs ) 
     573            ! 
     574            prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly 
     575            ! 
     576         END_2D 
     577         ! 
    615578         ! 
    616579      END SELECT 
    617580      ! 
    618       IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 
     581      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 
    619582      ! 
    620583      IF( ln_timing )   CALL timing_stop('eos2d') 
     
    648611      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    649612         ! 
    650          DO jk = 1, jpkm1 
    651             DO jj = 1, jpj 
    652                DO ji = 1, jpi 
    653                   ! 
    654                   zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
    655                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    656                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    657                   ztm = tmask(ji,jj,jk)                                         ! tmask 
    658                   ! 
    659                   ! alpha 
    660                   zn3 = ALP003 
    661                   ! 
    662                   zn2 = ALP012*zt + ALP102*zs+ALP002 
    663                   ! 
    664                   zn1 = ((ALP031*zt   & 
    665                      &   + ALP121*zs+ALP021)*zt   & 
    666                      &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
    667                      &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
    668                      ! 
    669                   zn0 = ((((ALP050*zt   & 
    670                      &   + ALP140*zs+ALP040)*zt   & 
    671                      &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
    672                      &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
    673                      &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
    674                      &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
    675                      ! 
    676                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    677                   ! 
    678                   pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 
    679                   ! 
    680                   ! beta 
    681                   zn3 = BET003 
    682                   ! 
    683                   zn2 = BET012*zt + BET102*zs+BET002 
    684                   ! 
    685                   zn1 = ((BET031*zt   & 
    686                      &   + BET121*zs+BET021)*zt   & 
    687                      &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
    688                      &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
    689                      ! 
    690                   zn0 = ((((BET050*zt   & 
    691                      &   + BET140*zs+BET040)*zt   & 
    692                      &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
    693                      &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
    694                      &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
    695                      &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
    696                      ! 
    697                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    698                   ! 
    699                   pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 
    700                   ! 
    701                END DO 
    702             END DO 
    703          END DO 
     613         DO_3D_11_11( 1, jpkm1 ) 
     614            ! 
     615            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
     616            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     617            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     618            ztm = tmask(ji,jj,jk)                                         ! tmask 
     619            ! 
     620            ! alpha 
     621            zn3 = ALP003 
     622            ! 
     623            zn2 = ALP012*zt + ALP102*zs+ALP002 
     624            ! 
     625            zn1 = ((ALP031*zt   & 
     626               &   + ALP121*zs+ALP021)*zt   & 
     627               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
     628               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
     629               ! 
     630            zn0 = ((((ALP050*zt   & 
     631               &   + ALP140*zs+ALP040)*zt   & 
     632               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
     633               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
     634               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
     635               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
     636               ! 
     637            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     638            ! 
     639            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 
     640            ! 
     641            ! beta 
     642            zn3 = BET003 
     643            ! 
     644            zn2 = BET012*zt + BET102*zs+BET002 
     645            ! 
     646            zn1 = ((BET031*zt   & 
     647               &   + BET121*zs+BET021)*zt   & 
     648               &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
     649               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
     650               ! 
     651            zn0 = ((((BET050*zt   & 
     652               &   + BET140*zs+BET040)*zt   & 
     653               &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
     654               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
     655               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
     656               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
     657               ! 
     658            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     659            ! 
     660            pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 
     661            ! 
     662         END_3D 
    704663         ! 
    705664      CASE( np_seos )                  !==  simplified EOS  ==! 
    706665         ! 
    707          DO jk = 1, jpkm1 
    708             DO jj = 1, jpj 
    709                DO ji = 1, jpi 
    710                   zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
    711                   zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
    712                   zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point 
    713                   ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask 
    714                   ! 
    715                   zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
    716                   pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha 
    717                   ! 
    718                   zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
    719                   pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta 
    720                   ! 
    721                END DO 
    722             END DO 
    723          END DO 
     666         DO_3D_11_11( 1, jpkm1 ) 
     667            zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
     668            zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
     669            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point 
     670            ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask 
     671            ! 
     672            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
     673            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha 
     674            ! 
     675            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
     676            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta 
     677            ! 
     678         END_3D 
    724679         ! 
    725680      CASE( np_leos )                  !==  linear ISOMIP EOS  ==! 
    726681         ! 
    727          DO jk = 1, jpkm1 
    728             DO jj = 1, jpj 
    729                DO ji = 1, jpi 
    730                   zt  = pts (ji,jj,jk,jp_tem) - (-1._wp) 
    731                   zs  = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
    732                   zh  = gdept(ji,jj,jk,Kmm)                 ! depth in meters at t-point 
    733                   ztm = tmask(ji,jj,jk)                   ! land/sea bottom mask = surf. mask 
    734                   ! 
    735                   zn  = rn_a0 * rho0 
    736                   pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha 
    737                   ! 
    738                   zn  = rn_b0 * rho0 
    739                   pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta 
    740                   ! 
    741                END DO 
    742             END DO 
    743          END DO 
     682         DO_3D_11_11( 1, jpkm1 ) 
     683            zt  = pts (ji,jj,jk,jp_tem) - (-1._wp) 
     684            zs  = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
     685            zh  = gdept(ji,jj,jk,Kmm)                 ! depth in meters at t-point 
     686            ztm = tmask(ji,jj,jk)                   ! land/sea bottom mask = surf. mask 
     687            ! 
     688            zn  = rn_a0 * rho0 
     689            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha 
     690            ! 
     691            zn  = rn_b0 * rho0 
     692            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta 
     693            ! 
     694         END_3D 
    744695         ! 
    745696      CASE DEFAULT 
     
    749700      END SELECT 
    750701      ! 
    751       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', & 
    752          &                       tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk ) 
     702      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', & 
     703         &                                  tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk ) 
    753704      ! 
    754705      IF( ln_timing )   CALL timing_stop('rab_3d') 
     
    783734      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    784735         ! 
    785          DO jj = 1, jpjm1 
    786             DO ji = 1, fs_jpim1   ! vector opt. 
    787                ! 
    788                zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
    789                zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
    790                zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    791                ! 
    792                ! alpha 
    793                zn3 = ALP003 
    794                ! 
    795                zn2 = ALP012*zt + ALP102*zs+ALP002 
    796                ! 
    797                zn1 = ((ALP031*zt   & 
    798                   &   + ALP121*zs+ALP021)*zt   & 
    799                   &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
    800                   &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
    801                   ! 
    802                zn0 = ((((ALP050*zt   & 
    803                   &   + ALP140*zs+ALP040)*zt   & 
    804                   &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
    805                   &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
    806                   &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
    807                   &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
    808                   ! 
    809                zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    810                ! 
    811                pab(ji,jj,jp_tem) = zn * r1_rho0 
    812                ! 
    813                ! beta 
    814                zn3 = BET003 
    815                ! 
    816                zn2 = BET012*zt + BET102*zs+BET002 
    817                ! 
    818                zn1 = ((BET031*zt   & 
    819                   &   + BET121*zs+BET021)*zt   & 
    820                   &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
    821                   &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
    822                   ! 
    823                zn0 = ((((BET050*zt   & 
    824                   &   + BET140*zs+BET040)*zt   & 
    825                   &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
    826                   &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
    827                   &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
    828                   &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
    829                   ! 
    830                zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    831                ! 
    832                pab(ji,jj,jp_sal) = zn / zs * r1_rho0 
    833                ! 
    834                ! 
    835             END DO 
    836          END DO 
    837          !                            ! Lateral boundary conditions 
    838          CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                     
     736         DO_2D_11_11 
     737            ! 
     738            zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
     739            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
     740            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     741            ! 
     742            ! alpha 
     743            zn3 = ALP003 
     744            ! 
     745            zn2 = ALP012*zt + ALP102*zs+ALP002 
     746            ! 
     747            zn1 = ((ALP031*zt   & 
     748               &   + ALP121*zs+ALP021)*zt   & 
     749               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
     750               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
     751               ! 
     752            zn0 = ((((ALP050*zt   & 
     753               &   + ALP140*zs+ALP040)*zt   & 
     754               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
     755               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
     756               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
     757               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
     758               ! 
     759            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     760            ! 
     761            pab(ji,jj,jp_tem) = zn * r1_rho0 
     762            ! 
     763            ! beta 
     764            zn3 = BET003 
     765            ! 
     766            zn2 = BET012*zt + BET102*zs+BET002 
     767            ! 
     768            zn1 = ((BET031*zt   & 
     769               &   + BET121*zs+BET021)*zt   & 
     770               &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
     771               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
     772               ! 
     773            zn0 = ((((BET050*zt   & 
     774               &   + BET140*zs+BET040)*zt   & 
     775               &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
     776               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
     777               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
     778               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
     779               ! 
     780            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     781            ! 
     782            pab(ji,jj,jp_sal) = zn / zs * r1_rho0 
     783            ! 
     784            ! 
     785         END_2D 
    839786         ! 
    840787      CASE( np_seos )                  !==  simplified EOS  ==! 
    841788         ! 
    842          DO jj = 1, jpjm1 
    843             DO ji = 1, fs_jpim1   ! vector opt. 
    844                ! 
    845                zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
    846                zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
    847                zh    = pdep (ji,jj)                   ! depth at the partial step level 
    848                ! 
    849                zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
    850                pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha 
    851                ! 
    852                zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
    853                pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta 
    854                ! 
    855             END DO 
    856          END DO 
    857          !                            ! Lateral boundary conditions 
    858          CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                     
     789         DO_2D_11_11 
     790            ! 
     791            zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
     792            zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
     793            zh    = pdep (ji,jj)                   ! depth at the partial step level 
     794            ! 
     795            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
     796            pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha 
     797            ! 
     798            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
     799            pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta 
     800            ! 
     801         END_2D 
    859802         ! 
    860803      CASE( np_leos )                  !==  linear ISOMIP EOS  ==! 
    861804         ! 
    862          DO jj = 1, jpjm1 
    863             DO ji = 1, fs_jpim1   ! vector opt. 
    864                ! 
    865                zt    = pts  (ji,jj,jp_tem) - (-1._wp)   ! pot. temperature anomaly (t-T0) 
    866                zs    = pts  (ji,jj,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
    867                zh    = pdep (ji,jj)                   ! depth at the partial step level 
    868                ! 
    869                zn  = rn_a0 * rho0 
    870                pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha 
    871                ! 
    872                zn  = rn_b0 * rho0 
    873                pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta 
    874                ! 
    875             END DO 
    876          END DO 
    877          ! 
    878          CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                    ! Lateral boundary conditions 
     805         DO_2D_11_11 
     806            ! 
     807            zt    = pts  (ji,jj,jp_tem) - (-1._wp)   ! pot. temperature anomaly (t-T0) 
     808            zs    = pts  (ji,jj,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
     809            zh    = pdep (ji,jj)                   ! depth at the partial step level 
     810            ! 
     811            zn  = rn_a0 * rho0 
     812            pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha 
     813            ! 
     814            zn  = rn_b0 * rho0 
     815            pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta 
     816            ! 
     817         END_2D 
    879818         ! 
    880819      CASE DEFAULT 
     
    884823      END SELECT 
    885824      ! 
    886       IF(ln_ctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', & 
    887          &                       tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' ) 
     825      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', & 
     826         &                                  tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' ) 
    888827      ! 
    889828      IF( ln_timing )   CALL timing_stop('rab_2d') 
     
    1026965      IF( ln_timing )   CALL timing_start('bn2') 
    1027966      ! 
    1028       DO jk = 2, jpkm1           ! interior points only (2=< jk =< jpkm1 ) 
    1029          DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90 
    1030             DO ji = 1, jpi 
    1031                zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   & 
    1032                   &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) )  
    1033                   ! 
    1034                zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw  
    1035                zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 
    1036                ! 
    1037                pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
    1038                   &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
    1039                   &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
    1040             END DO 
    1041          END DO 
    1042       END DO 
    1043       ! 
    1044       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', kdim=jpk ) 
     967      DO_3D_11_11( 2, jpkm1 ) 
     968         zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   & 
     969            &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) )  
     970            ! 
     971         zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw  
     972         zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 
     973         ! 
     974         pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
     975            &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
     976            &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
     977      END_3D 
     978      ! 
     979      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', kdim=jpk ) 
    1045980      ! 
    1046981      IF( ln_timing )   CALL timing_stop('bn2') 
     
    10781013      z1_T0   = 1._wp/40._wp 
    10791014      ! 
    1080       DO jj = 1, jpj 
    1081          DO ji = 1, jpi 
    1082             ! 
    1083             zt  = ctmp   (ji,jj) * z1_T0 
    1084             zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 ) 
    1085             ztm = tmask(ji,jj,1) 
    1086             ! 
    1087             zn = ((((-2.1385727895e-01_wp*zt   & 
    1088                &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   & 
    1089                &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   & 
    1090                &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   & 
    1091                &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   & 
    1092                &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   & 
    1093                &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   & 
    1094                &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp 
    1095                ! 
    1096             zd = (2.0035003456_wp*zt   & 
    1097                &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   & 
    1098                &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp 
    1099                ! 
    1100             ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm 
    1101                ! 
    1102          END DO 
    1103       END DO 
     1015      DO_2D_11_11 
     1016         ! 
     1017         zt  = ctmp   (ji,jj) * z1_T0 
     1018         zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 ) 
     1019         ztm = tmask(ji,jj,1) 
     1020         ! 
     1021         zn = ((((-2.1385727895e-01_wp*zt   & 
     1022            &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   & 
     1023            &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   & 
     1024            &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   & 
     1025            &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   & 
     1026            &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   & 
     1027            &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   & 
     1028            &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp 
     1029            ! 
     1030         zd = (2.0035003456_wp*zt   & 
     1031            &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   & 
     1032            &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp 
     1033            ! 
     1034         ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm 
     1035            ! 
     1036      END_2D 
    11041037      ! 
    11051038      IF( ln_timing )   CALL timing_stop('eos_pt_from_ct') 
     
    11331066         ! 
    11341067         z1_S0 = 1._wp / 35.16504_wp 
    1135          DO jj = 1, jpj 
    1136             DO ji = 1, jpi 
    1137                zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity 
    1138                ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 
    1139                   &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 
    1140             END DO 
    1141          END DO 
     1068         DO_2D_11_11 
     1069            zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity 
     1070            ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 
     1071               &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 
     1072         END_2D 
    11421073         ptf(:,:) = ptf(:,:) * psal(:,:) 
    11431074         ! 
    11441075         IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:) 
    11451076         ! 
    1146       CASE ( np_eos80, np_leos )                !==  PT,SP (UNESCO formulation)  ==! 
     1077      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==! 
    11471078         ! 
    11481079         ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   & 
     
    11901121         IF( PRESENT( pdep ) )   ptf = ptf - 7.53e-4 * pdep 
    11911122         ! 
    1192       CASE ( np_eos80, np_leos )                !==  PT,SP (UNESCO formulation)  ==! 
     1123      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==! 
    11931124         ! 
    11941125         ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal )   & 
     
    12421173      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    12431174         ! 
    1244          DO jk = 1, jpkm1 
    1245             DO jj = 1, jpj 
    1246                DO ji = 1, jpi 
    1247                   ! 
    1248                   zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
    1249                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    1250                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    1251                   ztm = tmask(ji,jj,jk)                                         ! tmask 
    1252                   ! 
    1253                   ! potential energy non-linear anomaly 
    1254                   zn2 = (PEN012)*zt   & 
    1255                      &   + PEN102*zs+PEN002 
    1256                      ! 
    1257                   zn1 = ((PEN021)*zt   & 
    1258                      &   + PEN111*zs+PEN011)*zt   & 
    1259                      &   + (PEN201*zs+PEN101)*zs+PEN001 
    1260                      ! 
    1261                   zn0 = ((((PEN040)*zt   & 
    1262                      &   + PEN130*zs+PEN030)*zt   & 
    1263                      &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   & 
    1264                      &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   & 
    1265                      &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000 
    1266                      ! 
    1267                   zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    1268                   ! 
    1269                   ppen(ji,jj,jk)  = zn * zh * r1_rho0 * ztm 
    1270                   ! 
    1271                   ! alphaPE non-linear anomaly 
    1272                   zn2 = APE002 
    1273                   ! 
    1274                   zn1 = (APE011)*zt   & 
    1275                      &   + APE101*zs+APE001 
    1276                      ! 
    1277                   zn0 = (((APE030)*zt   & 
    1278                      &   + APE120*zs+APE020)*zt   & 
    1279                      &   + (APE210*zs+APE110)*zs+APE010)*zt   & 
    1280                      &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000 
    1281                      ! 
    1282                   zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    1283                   !                               
    1284                   pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm 
    1285                   ! 
    1286                   ! betaPE non-linear anomaly 
    1287                   zn2 = BPE002 
    1288                   ! 
    1289                   zn1 = (BPE011)*zt   & 
    1290                      &   + BPE101*zs+BPE001 
    1291                      ! 
    1292                   zn0 = (((BPE030)*zt   & 
    1293                      &   + BPE120*zs+BPE020)*zt   & 
    1294                      &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   & 
    1295                      &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000 
    1296                      ! 
    1297                   zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    1298                   !                               
    1299                   pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm 
    1300                   ! 
    1301                END DO 
    1302             END DO 
    1303          END DO 
     1175         DO_3D_11_11( 1, jpkm1 ) 
     1176            ! 
     1177            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
     1178            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     1179            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     1180            ztm = tmask(ji,jj,jk)                                         ! tmask 
     1181            ! 
     1182            ! potential energy non-linear anomaly 
     1183            zn2 = (PEN012)*zt   & 
     1184               &   + PEN102*zs+PEN002 
     1185               ! 
     1186            zn1 = ((PEN021)*zt   & 
     1187               &   + PEN111*zs+PEN011)*zt   & 
     1188               &   + (PEN201*zs+PEN101)*zs+PEN001 
     1189               ! 
     1190            zn0 = ((((PEN040)*zt   & 
     1191               &   + PEN130*zs+PEN030)*zt   & 
     1192               &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   & 
     1193               &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   & 
     1194               &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000 
     1195               ! 
     1196            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     1197            ! 
     1198            ppen(ji,jj,jk)  = zn * zh * r1_rho0 * ztm 
     1199            ! 
     1200            ! alphaPE non-linear anomaly 
     1201            zn2 = APE002 
     1202            ! 
     1203            zn1 = (APE011)*zt   & 
     1204               &   + APE101*zs+APE001 
     1205               ! 
     1206            zn0 = (((APE030)*zt   & 
     1207               &   + APE120*zs+APE020)*zt   & 
     1208               &   + (APE210*zs+APE110)*zs+APE010)*zt   & 
     1209               &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000 
     1210               ! 
     1211            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     1212            !                               
     1213            pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm 
     1214            ! 
     1215            ! betaPE non-linear anomaly 
     1216            zn2 = BPE002 
     1217            ! 
     1218            zn1 = (BPE011)*zt   & 
     1219               &   + BPE101*zs+BPE001 
     1220               ! 
     1221            zn0 = (((BPE030)*zt   & 
     1222               &   + BPE120*zs+BPE020)*zt   & 
     1223               &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   & 
     1224               &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000 
     1225               ! 
     1226            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     1227            !                               
     1228            pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm 
     1229            ! 
     1230         END_3D 
    13041231         ! 
    13051232      CASE( np_seos )                !==  Vallis (2006) simplified EOS  ==! 
    13061233         ! 
    1307          DO jk = 1, jpkm1 
    1308             DO jj = 1, jpj 
    1309                DO ji = 1, jpi 
    1310                   zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0) 
    1311                   zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0) 
    1312                   zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point 
    1313                   ztm = tmask(ji,jj,jk)                ! tmask 
    1314                   zn  = 0.5_wp * zh * r1_rho0 * ztm 
    1315                   !                                    ! Potential Energy 
    1316                   ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn 
    1317                   !                                    ! alphaPE 
    1318                   pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn 
    1319                   pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn 
    1320                   ! 
    1321                END DO 
    1322             END DO 
    1323          END DO 
     1234         DO_3D_11_11( 1, jpkm1 ) 
     1235            zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0) 
     1236            zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0) 
     1237            zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point 
     1238            ztm = tmask(ji,jj,jk)                ! tmask 
     1239            zn  = 0.5_wp * zh * r1_rho0 * ztm 
     1240            !                                    ! Potential Energy 
     1241            ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn 
     1242            !                                    ! alphaPE 
     1243            pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn 
     1244            pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn 
     1245            ! 
     1246         END_3D 
    13241247         ! 
    13251248      CASE( np_leos )                !==  linear ISOMIP EOS  ==! 
    13261249         ! 
    1327          DO jk = 1, jpkm1 
    1328             DO jj = 1, jpj 
    1329                DO ji = 1, jpi 
    1330                   zt  = pts(ji,jj,jk,jp_tem) - (-1._wp)  ! temperature anomaly (t-T0) 
    1331                   zs = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
    1332                   zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters  at t-point 
    1333                   ztm = tmask(ji,jj,jk)                  ! tmask 
    1334                   zn  = 0.5_wp * zh * r1_rho0 * ztm 
    1335                   !                                    ! Potential Energy 
    1336                   ppen(ji,jj,jk) = 0. 
    1337                   !                                    ! alphaPE 
    1338                   pab_pe(ji,jj,jk,jp_tem) = 0. 
    1339                   pab_pe(ji,jj,jk,jp_sal) = 0. 
    1340                   ! 
    1341                END DO 
    1342             END DO 
    1343          END DO 
     1250         DO_3D_11_11( 1, jpkm1 ) 
     1251            zt  = pts(ji,jj,jk,jp_tem) - (-1._wp)  ! temperature anomaly (t-T0) 
     1252            zs = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0) 
     1253            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters  at t-point 
     1254            ztm = tmask(ji,jj,jk)                  ! tmask 
     1255            zn  = 0.5_wp * zh * r1_rho0 * ztm 
     1256            !                                    ! Potential Energy 
     1257            ppen(ji,jj,jk) = 0. 
     1258            !                                    ! alphaPE 
     1259            pab_pe(ji,jj,jk,jp_tem) = 0. 
     1260            pab_pe(ji,jj,jk,jp_sal) = 0. 
     1261            ! 
     1262         END_3D 
    13441263         ! 
    13451264      CASE DEFAULT 
     
    13651284      INTEGER  ::   ioptio   ! local integer 
    13661285      !! 
    1367       NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS   , ln_LEOS, & 
    1368          &             rn_a0    , rn_b0   , rn_lambda1, rn_mu1 , & 
    1369          &                                  rn_lambda2, rn_mu2 , rn_nu 
    1370       !!---------------------------------------------------------------------- 
    1371       ! 
    1372       REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state 
     1286      NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, ln_LEOS, rn_a0, rn_b0, & 
     1287         &             rn_lambda1, rn_mu1, rn_lambda2, rn_mu2, rn_nu 
     1288      !!---------------------------------------------------------------------- 
     1289      ! 
    13731290      READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 ) 
    13741291901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist' ) 
    13751292      ! 
    1376       REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state 
    13771293      READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 ) 
    13781294902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist' ) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP+/MY_SRC/isfcavgam.F90

    r12077 r13159  
    9191         pgs(:,:) = rn_gammas0 
    9292      CASE ( 'vel' ) ! gamma is proportional to u* 
    93          CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, r_ke0_top,               pgt, pgs ) 
     93         CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, rn_vtide**2,               pgt, pgs ) 
    9494      CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u* 
    95          CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pgt, pgs ) 
     95         CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, rn_vtide**2, pqoce, pqfwf, pgt, pgs ) 
    9696      CASE DEFAULT 
    9797         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)') 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP+/MY_SRC/isfstp.F90

    r12077 r13159  
    250250      IF ( l_isfoasis .AND. ln_isf ) THEN 
    251251         ! 
    252          CALL ctl_stop( ' ln_ctl and ice shelf not tested' ) 
     252         CALL ctl_stop( 'namelist combination ln_cpl and ln_isf not tested' ) 
    253253         ! 
    254254         ! NEMO coupled to ATMO model with isf cavity need oasis method for melt computation  
     
    291291      !!---------------------------------------------------------------------- 
    292292      ! 
    293       REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs  
    294293      READ  ( numnam_ref, namisf, IOSTAT = ios, ERR = 901) 
    295294901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist' ) 
    296295      ! 
    297       REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs 
    298296      READ  ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 ) 
    299297902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist' ) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP+/MY_SRC/istate.F90

    r12353 r13159  
    4141   PUBLIC   istate_init   ! routine called by step.F90 
    4242 
     43   !! * Substitutions 
     44#  include "do_loop_substitute.h90" 
    4345   !!---------------------------------------------------------------------- 
    4446   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7577      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
    7678      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk 
    77       ts   (:,:,:,:,Kaa) = 0._wp                               ! set one for all to 0 at level jpk 
     79      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk 
    7880      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk 
    7981#if defined key_agrif 
     
    9092         !                                    ! --------------- 
    9193         numror = 0                           ! define numror = 0 -> no restart file to read 
    92          neuler = 0                           ! Set time-step indicator at nit000 (euler forward) 
     94         l_1st_euler = .true.                 ! Set time-step indicator at nit000 (euler forward) 
    9395         CALL day_init                        ! model calendar (using both namelist and restart infos) 
    9496         !                                    ! Initialization of ocean to zero 
     
    103105               ! Apply minimum wetdepth criterion 
    104106               ! 
    105                DO jj = 1,jpj 
    106                   DO ji = 1,jpi 
    107                      IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN 
    108                         ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 
    109                      ENDIF 
    110                   END DO 
    111                END DO  
     107               DO_2D_11_11 
     108                  IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN 
     109                     ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 
     110                  ENDIF 
     111               END_2D 
    112112            ENDIF  
    113113            uu  (:,:,:,Kbb) = 0._wp 
     
    159159      ! 
    160160!!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked 
    161       DO jk = 1, jpkm1 
    162          DO jj = 1, jpj 
    163             DO ji = 1, jpi 
    164                uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk) 
    165                vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk) 
    166                ! 
    167                uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk) 
    168                vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk) 
    169             END DO 
    170          END DO 
    171       END DO 
     161      DO_3D_11_11( 1, jpkm1 ) 
     162         uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk) 
     163         vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk) 
     164         ! 
     165         uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk) 
     166         vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk) 
     167      END_3D 
    172168      ! 
    173169      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP+/MY_SRC/sbcfwb.F90

    r12489 r13159  
    151151         ENDIF    
    152152         !                                         ! Update fwfold if new year start 
    153          ikty = 365 * 86400 / rn_Dt               !!bug  use of 365 days leap year or 360d year !!!!!!! 
     153         ikty = 365 * 86400 / rn_Dt                  !!bug  use of 365 days leap year or 360d year !!!!!!! 
    154154         IF( MOD( kt, ikty ) == 0 ) THEN 
    155155            a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP+/MY_SRC/tradmp.F90

    r12353 r13159  
    5151   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   resto    !: restoring coeff. on T and S (s-1) 
    5252 
     53   !! * Substitutions 
     54#  include "do_loop_substitute.h90" 
    5355   !!---------------------------------------------------------------------- 
    5456   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    110112      CASE( 0 )                        !*  newtonian damping throughout the water column  *! 
    111113         DO jn = 1, jpts 
    112             DO jk = 1, jpkm1 
    113                DO jj = 2, jpjm1 
    114                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    115                      pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           & 
    116                         &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) ) 
    117                   END DO 
    118                END DO 
    119             END DO 
     114            DO_3D_00_00( 1, jpkm1 ) 
     115               pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           & 
     116                  &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) ) 
     117            END_3D 
    120118         END DO 
    121119         ! 
    122120      CASE ( 1 )                       !*  no damping in the turbocline (avt > 5 cm2/s)  *! 
    123          DO jk = 1, jpkm1 
    124             DO jj = 2, jpjm1 
    125                DO ji = fs_2, fs_jpim1   ! vector opt. 
    126                   IF( avt(ji,jj,jk) <= avt_c ) THEN 
    127                      pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
    128                         &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) ) 
    129                      pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   & 
    130                         &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) ) 
    131                   ENDIF 
    132                END DO 
    133             END DO 
    134          END DO 
     121         DO_3D_00_00( 1, jpkm1 ) 
     122            IF( avt(ji,jj,jk) <= avt_c ) THEN 
     123               pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
     124                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) ) 
     125               pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   & 
     126                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) ) 
     127            ENDIF 
     128         END_3D 
    135129         ! 
    136130      CASE ( 2 )                       !*  no damping in the mixed layer   *! 
    137          DO jk = 1, jpkm1 
    138             DO jj = 2, jpjm1 
    139                DO ji = fs_2, fs_jpim1   ! vector opt. 
    140                   IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN 
    141                      pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
    142                         &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) ) 
    143                      pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   & 
    144                         &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) ) 
    145                   ENDIF 
    146                END DO 
    147             END DO 
    148          END DO 
     131         DO_3D_00_00( 1, jpkm1 ) 
     132            IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN 
     133               pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
     134                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) ) 
     135               pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   & 
     136                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) ) 
     137            ENDIF 
     138         END_3D 
    149139         ! 
    150140      END SELECT 
     
    157147      ENDIF 
    158148      !                           ! Control print 
    159       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   & 
    160          &                       tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     149      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   & 
     150         &                                  tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    161151      ! 
    162152      IF( ln_timing )   CALL timing_stop('tra_dmp') 
     
    178168      !!---------------------------------------------------------------------- 
    179169      ! 
    180       REWIND( numnam_ref )   ! Namelist namtra_dmp in reference namelist : T & S relaxation 
    181170      READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901) 
    182171901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_dmp in reference namelist' ) 
    183172      ! 
    184       REWIND( numnam_cfg )   ! Namelist namtra_dmp in configuration namelist : T & S relaxation 
    185173      READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 ) 
    186174902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist' ) 
Note: See TracChangeset for help on using the changeset viewer.