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 13214 for NEMO/trunk – NEMO

Changeset 13214 for NEMO/trunk


Ignore:
Timestamp:
2020-07-02T11:09:01+02:00 (4 years ago)
Author:
smasson
Message:

trunk: Mid-year merge, merge back dev_r12563_ASINTER-06_ABL_improvement

Location:
NEMO/trunk
Files:
13 edited
1 copied

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-oce.xml

    r12063 r13214  
    5656     <field field_ref="t_abl" /> 
    5757     <field field_ref="q_abl" /> 
     58     <field field_ref="uvz1_abl" /> 
     59     <field field_ref="tz1_abl" /> 
     60     <field field_ref="qz1_abl" /> 
     61     <field field_ref="uvz1_dta" /> 
     62     <field field_ref="tz1_dta" /> 
     63     <field field_ref="qz1_dta" /> 
    5864     <field field_ref="pblh" /> 
    5965   </file> 
  • NEMO/trunk/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg

    r12994 r13214  
    121121/ 
    122122!----------------------------------------------------------------------- 
     123&namsbc_abl    !   Atmospheric Boundary Layer formulation           (ln_abl = T) 
     124!----------------------------------------------------------------------- 
     125/ 
     126!----------------------------------------------------------------------- 
    123127&namtra_qsr    !   penetrative solar radiation                          (ln_traqsr =T) 
    124128!----------------------------------------------------------------------- 
  • NEMO/trunk/cfgs/SHARED/field_def_nemo-oce.xml

    r13132 r13214  
    455455          <field id="t_dta"      long_name="DTA potential temperature"     standard_name="dta_theta"      unit="K"        /> 
    456456          <field id="q_dta"      long_name="DTA specific humidity"         standard_name="dta_qspe"       unit="kg/kg"    /> 
    457           <field id="coeft"      long_name="ABL nudging coefficient"       standard_name="coeft"          unit=""         /> 
     457          <field id="u_geo"      long_name="GEO i-horizontal velocity"     standard_name="geo_x_velocity" unit="m/s"      /> 
     458          <field id="v_geo"      long_name="GEO j-horizontal velocity"     standard_name="geo_y_velocity" unit="m/s"      /> 
    458459          <field id="tke_abl"    long_name="ABL turbulent kinetic energy"  standard_name="abl_tke"        unit="m2/s2"    /> 
    459460          <field id="avm_abl"    long_name="ABL turbulent viscosity"       standard_name="abl_avm"        unit="m2/s"     /> 
    460461          <field id="avt_abl"    long_name="ABL turbulent diffusivity"     standard_name="abl_avt"        unit="m2/s"     /> 
    461           <field id="mxl_abl"    long_name="ABL mixing length"             standard_name="abl_mxl"        unit="m"        /> 
     462          <field id="mxlm_abl"   long_name="ABL master mixing length"      standard_name="abl_mxlm"       unit="m"        /> 
     463          <field id="mxld_abl"   long_name="ABL dissipative mixing length" standard_name="abl_mxld"       unit="m"        /> 
    462464   </field_group> 
    463465 
  • NEMO/trunk/cfgs/SHARED/namelist_ref

    r13208 r13214  
    279279      ln_humi_dpt = .false. !  humidity "sn_humi" is dew-point temperature [K] 
    280280      ln_humi_rlh = .false. !  humidity "sn_humi" is relative humidity     [%] 
     281      ln_tpot     = .true.  !!GS: compute potential temperature or not 
    281282   ! 
    282283   cn_dir      = './'      !  root directory for the bulk data location 
     
    309310   cn_ablrst_outdir = "."             !  directory to write output abl restarts 
    310311 
     312   ln_rstart_abl  = .false. 
    311313   ln_hpgls_frc   = .false. 
    312314   ln_geos_winds  = .false. 
    313    nn_dyn_restore = 2         ! restoring option for dynamical ABL variables: = 0 no restoring 
     315   ln_smth_pblh   = .false. 
     316   nn_dyn_restore = 0         ! restoring option for dynamical ABL variables: = 0 no restoring 
    314317                              !                                               = 1 equatorial restoring 
    315318                              !                                               = 2 global restoring 
    316    rn_ldyn_min   =  4.5       !  magnitude of the nudging on ABL dynamics at the bottom of the ABL   [hour] 
    317    rn_ldyn_max   =  1.5       !  magnitude of the nudging on ABL dynamics at the top of the ABL   [hour] 
    318    rn_ltra_min   =  4.5       !  magnitude of the nudging on ABL tracers  at the bottom of the ABL   [hour] 
    319    rn_ltra_max   =  1.5       !  magnitude of the nudging on ABL tracers  at the top of the ABL   [hour] 
     319   rn_ldyn_min   =  4.5       ! dynamics nudging magnitude inside the ABL [hour] (~3 rn_Dt) 
     320   rn_ldyn_max   =  1.5       ! dynamics nudging magnitude above  the ABL [hour] (~1 rn_Dt) 
     321   rn_ltra_min   =  4.5       ! tracers  nudging magnitude inside the ABL [hour] (~3 rn_Dt) 
     322   rn_ltra_max   =  1.5       ! tracers  nudging magnitude above  the ABL [hour] (~1 rn_Dt) 
    320323   nn_amxl       =  0         ! mixing length: = 0 Deardorff 80 length-scale 
    321324                              !                = 1 length-scale based on the distance to the PBL height 
    322325                              !                = 2 Bougeault & Lacarrere 89 length-scale 
    323    rn_Cm         = 0.0667     ! 0.126 in MesoNH 
    324    rn_Ct         = 0.1667     ! 0.143 in MesoNH 
    325    rn_Ce         = 0.4        ! 0.4   in MesoNH 
    326    rn_Ceps       = 0.7        ! 0.85  in MesoNH 
    327    rn_Rod        = 0.15       ! c0 in RMCA17 mixing length formulation (not yet implemented) 
    328    rn_Ric        = 0.139      !  Critical Richardson number (to compute PBL height and diffusivities) 
     326                              ! CBR00  ! CCH02  ! MesoNH ! 
     327   rn_Cm          = 0.0667    ! 0.0667 ! 0.1260 ! 0.1260 ! 
     328   rn_Ct          = 0.1667    ! 0.1667 ! 0.1430 ! 0.1430 ! 
     329   rn_Ce          = 0.40      ! 0.40   ! 0.34   ! 0.40   ! 
     330   rn_Ceps        = 0.700     ! 0.700  ! 0.845  ! 0.850  ! 
     331   rn_Ric         = 0.139     ! 0.139  ! 0.143  !   ?    ! Critical Richardson number (to compute PBL height and diffusivities) 
     332   rn_Rod         = 0.15      ! c0 in RMCA17 mixing length formulation (not yet implemented) 
    329333/ 
    330334!----------------------------------------------------------------------- 
  • NEMO/trunk/src/ABL/abl.F90

    r12489 r13214  
    2929   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   avm_abl      !: turbulent viscosity   [m2/s] 
    3030   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   avt_abl      !: turbulent diffusivity [m2/s] 
    31    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   mxl_abl      !: mixing length         [m] 
     31   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   mxld_abl     !: dissipative mixing length    [m] 
     32   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   mxlm_abl     !: master mixing length         [m] 
    3233   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:)   ::   tke_abl      !: turbulent kinetic energy [m2/s2] 
    3334   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   fft_abl      !: Coriolis parameter    [1/s] 
     
    5556      !!---------------------------------------------------------------------- 
    5657      ! 
    57       ALLOCATE( u_abl  (1:jpi,1:jpj,1:jpka,jptime), & 
    58          &      v_abl  (1:jpi,1:jpj,1:jpka,jptime), & 
    59          &      tq_abl (1:jpi,1:jpj,1:jpka,jptime,jptq), & 
    60          &      avm_abl(1:jpi,1:jpj,1:jpka), & 
    61          &      avt_abl(1:jpi,1:jpj,1:jpka), & 
    62          &      mxl_abl(1:jpi,1:jpj,1:jpka), &          
    63          &      tke_abl(1:jpi,1:jpj,1:jpka,jptime), & 
    64          &      fft_abl(1:jpi,1:jpj), & 
    65          &      pblh   (1:jpi,1:jpj), &          
    66          &      msk_abl(1:jpi,1:jpj), & 
    67          &      rest_eq(1:jpi,1:jpj), &          
    68          &      e3t_abl(1:jpka), e3w_abl(1:jpka), ght_abl(1:jpka), ghw_abl(1:jpka),                 STAT=ierr ) 
     58      ALLOCATE( u_abl   (1:jpi,1:jpj,1:jpka,jptime     ), & 
     59         &      v_abl   (1:jpi,1:jpj,1:jpka,jptime     ), & 
     60         &      tq_abl  (1:jpi,1:jpj,1:jpka,jptime,jptq), & 
     61         &      tke_abl (1:jpi,1:jpj,1:jpka,jptime     ), & 
     62         &      avm_abl (1:jpi,1:jpj,1:jpka            ), & 
     63         &      avt_abl (1:jpi,1:jpj,1:jpka            ), & 
     64         &      mxld_abl(1:jpi,1:jpj,1:jpka            ), & 
     65         &      mxlm_abl(1:jpi,1:jpj,1:jpka            ), & 
     66         &      fft_abl (1:jpi,1:jpj                   ), & 
     67         &      pblh    (1:jpi,1:jpj                   ), & 
     68         &      msk_abl (1:jpi,1:jpj                   ), & 
     69         &      rest_eq (1:jpi,1:jpj                   ), & 
     70         &      e3t_abl (1:jpka), e3w_abl(1:jpka)       , & 
     71         &      ght_abl (1:jpka), ghw_abl(1:jpka)       , STAT=ierr ) 
    6972         ! 
    7073      abl_alloc = ierr 
  • NEMO/trunk/src/ABL/ablmod.F90

    r13208 r13214  
    22   !!====================================================================== 
    33   !!                       ***  MODULE  ablmod  *** 
    4    !! Surface module :  ABL computation to provide atmospheric data  
     4   !! Surface module :  ABL computation to provide atmospheric data 
    55   !!                   for surface fluxes computation 
    66   !!====================================================================== 
    77   !! History :  3.6  ! 2019-03  (F. Lemarié & G. Samson)  Original code 
    88   !!---------------------------------------------------------------------- 
    9     
     9 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   abl_stp       : ABL single column model 
     
    1616 
    1717   USE phycst         ! physical constants 
    18    USE dom_oce, ONLY  : tmask   
     18   USE dom_oce, ONLY  : tmask 
    1919   USE sbc_oce, ONLY  : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1, rhoa 
    20    USE sbcblk, ONLY   : rn_efac 
     20   USE sbcblk         ! use rn_efac, cdn_oce 
    2121   USE sbcblk_phy     ! use some physical constants for flux computation 
    2222   ! 
     
    3030 
    3131   PUBLIC   abl_stp   ! called by sbcabl.F90 
    32    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ustar2 
     32   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ustar2, zrough 
    3333   !! * Substitutions 
    3434#  include "do_loop_substitute.h90" 
     
    3838 
    3939!=================================================================================================== 
    40    SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, &                            ! in 
    41               &            pu_dta, pv_dta, pt_dta, pq_dta,    &  
     40   SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq,            &     ! in 
     41              &            pu_dta, pv_dta, pt_dta, pq_dta,    & 
    4242              &            pslp_dta, pgu_dta, pgv_dta,        & 
    43               &            pcd_du, psen, pevp,                &                            ! in/out 
    44               &            pwndm, ptaui, ptauj, ptaum         &  
    45 #if defined key_si3           
    46               &          , ptm_su,pssu_ice,pssv_ice,pssq_ice,pcd_du_ice  &  
    47               &          , psen_ice, pevp_ice, pwndm_ice, pfrac_oce      &  
    48               &          , ptaui_ice, ptauj_ice                          & 
    49 #endif            
    50            &   ) 
     43              &            pcd_du, psen, pevp,                &     ! in/out 
     44              &            pwndm, ptaui, ptauj, ptaum         & 
     45#if defined key_si3 
     46              &          , ptm_su, pssu_ice, pssv_ice         & 
     47              &          , pssq_ice, pcd_du_ice, psen_ice     & 
     48              &          , pevp_ice, pwndm_ice, pfrac_oce     & 
     49              &          , ptaui_ice, ptauj_ice               & 
     50#endif 
     51              &      ) 
    5152!--------------------------------------------------------------------------------------------------- 
    5253 
     
    5455      !!                    ***  ROUTINE abl_stp *** 
    5556      !! 
    56       !! ** Purpose :   Time-integration of the ABL model  
     57      !! ** Purpose :   Time-integration of the ABL model 
    5758      !! 
    58       !! ** Method  :   Compute atmospheric variables : vertical turbulence  
     59      !! ** Method  :   Compute atmospheric variables : vertical turbulence 
    5960      !!                             + Coriolis term + newtonian relaxation 
    60       !!                 
     61      !! 
    6162      !! ** Action  : - Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh 
    6263      !!              - Advance tracers to time n+1 (Euler backward scheme) 
     
    7071      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psst       ! sea-surface temperature [Celsius] 
    7172      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu       ! sea-surface u (U-point) 
    72       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv       ! sea-surface v (V-point)  
     73      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv       ! sea-surface v (V-point) 
    7374      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq       ! sea-surface humidity 
    7475      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pu_dta     ! large-scale windi 
     
    8283      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   psen       ! Ch x Du 
    8384      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pevp       ! Ce x Du 
    84       REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pwndm      ! ||uwnd||       
     85      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pwndm      ! ||uwnd|| 
    8586      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui      ! taux 
    8687      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj      ! tauy 
    87       REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaum      ! ||tau||  
    88       ! 
    89 #if defined key_si3     
     88      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaum      ! ||tau|| 
     89      ! 
     90#if defined key_si3 
    9091      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptm_su       ! ice-surface temperature [K] 
    9192      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu_ice     ! ice-surface u (U-point) 
    92       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv_ice     ! ice-surface v (V-point)       
    93       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq_ice     ! ice-surface humidity  
     93      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv_ice     ! ice-surface v (V-point) 
     94      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq_ice     ! ice-surface humidity 
    9495      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du_ice   ! Cd x Du over ice (T-point) 
    9596      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psen_ice     ! Ch x Du over ice (T-point) 
    9697      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pevp_ice     ! Ce x Du over ice (T-point) 
    97       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndm_ice    ! ||uwnd - uice||   
    98       !REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pfrac_oce     !!GS: out useless ? 
    99       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pfrac_oce      ! 
     98      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndm_ice    ! ||uwnd - uice|| 
     99      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pfrac_oce    ! ocean fraction 
    100100      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui_ice    ! ice-surface taux stress (U-point) 
    101       REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (V-point)      
    102 #endif      
    103       ! 
    104       REAL(wp), DIMENSION(1:jpi,1:jpj   )        ::   zwnd_i, zwnd_j 
    105       REAL(wp), DIMENSION(1:jpi,2:jpka  )        ::   zCF     
    106       REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka)    ::   z_cft      !--FL--to be removed after the test phase    
    107       ! 
    108       REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_a 
    109       REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_b 
    110       REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_c 
     101      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (V-point) 
     102#endif 
     103      ! 
     104      REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   zwnd_i, zwnd_j 
     105      REAL(wp), DIMENSION(1:jpi      ,2:jpka)    ::   zCF 
     106      ! 
     107      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_a 
     108      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_b 
     109      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_c 
    111110      ! 
    112111      INTEGER             ::   ji, jj, jk, jtra, jbak               ! dummy loop indices 
    113112      REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2 
    114113      REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce 
    115       ! 
    116       !!---------------------------------------------------------------------       
     114      LOGICAL             ::   SemiImp_Cor = .TRUE. 
     115      ! 
     116      !!--------------------------------------------------------------------- 
    117117      ! 
    118118      IF(lwp .AND. kt == nit000) THEN                  ! control print 
     
    120120         WRITE(numout,*) 'abl_stp : ABL time stepping' 
    121121         WRITE(numout,*) '~~~~~~' 
    122       ENDIF  
     122      ENDIF 
    123123      ! 
    124124      IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) ) 
    125       !! Compute ustar squared as Cd || Uatm-Uoce ||^2  
    126       !! needed for surface boundary condition of TKE  
     125      IF( kt == nit000 ) ALLOCATE ( zrough( 1:jpi, 1:jpj ) ) 
     126      !! Compute ustar squared as Cd || Uatm-Uoce ||^2 
     127      !! needed for surface boundary condition of TKE 
    127128      !! pwndm contains | U10m - U_oce | (see blk_oce_1 in sbcblk) 
    128129      DO_2D_11_11 
    129130         zzoce         = pCd_du    (ji,jj) * pwndm    (ji,jj) 
    130131#if defined key_si3 
    131          zzice         = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj)  
    132          ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice  
     132         zzice         = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj) 
     133         ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice 
    133134#else 
    134          ustar2(ji,jj) = zzoce    
     135         ustar2(ji,jj) = zzoce 
    135136#endif 
     137         zrough(ji,jj) = ght_abl(2) * EXP( - vkarmn / SQRT( MAX( Cdn_oce(ji,jj), 1.e-4 ) ) ) !<-- recover the value of z0 from Cdn_oce 
    136138      END_2D 
    137139      ! 
     
    140142      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    141143 
    142       CALL abl_zdf_tke( )                       !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj)  
    143           
     144      CALL abl_zdf_tke( )                       !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj) 
     145 
    144146      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    145147      !                            !  2 *** Advance tracers to time n+1 
    146148      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    147        
     149 
    148150      !------------- 
    149151      DO jj = 1, jpj    ! outer loop             !--> tq_abl computed on (1:jpi) x (1:jpj) 
    150       !-------------     
    151          ! Compute matrix elements for interior points  
     152      !------------- 
     153         ! Compute matrix elements for interior points 
    152154         DO jk = 3, jpkam1 
    153155            DO ji = 1, jpi   ! vector opt. 
    154                z_elem_a( ji,     jk              ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal 
    155                z_elem_c( ji,     jk              ) = - rDt_abl * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal        
    156                z_elem_b( ji,     jk              ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   !       diagonal            
    157             END DO 
    158          END DO   
    159          ! Boundary conditions   
    160          DO ji = 1, jpi   ! vector opt.  
    161             ! Neumann at the bottom           
    162             z_elem_a( ji,     2              ) = 0._wp 
    163             z_elem_c( ji,     2              ) = - rDt_abl * Avt_abl( ji, jj, 2   ) / e3w_abl( 2   )                                          
     156               z_elem_a( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal 
     157               z_elem_c( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal 
     158               z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   !       diagonal 
     159            END DO 
     160         END DO 
     161         ! Boundary conditions 
     162         DO ji = 1, jpi   ! vector opt. 
     163            ! Neumann at the bottom 
     164            z_elem_a( ji,    2 ) = 0._wp 
     165            z_elem_c( ji,    2 ) = - rDt_abl * Avt_abl( ji, jj,    2 ) / e3w_abl(    2 ) 
    164166            ! Homogeneous Neumann at the top 
    165             z_elem_a( ji,     jpka           ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )  
    166             z_elem_c( ji,     jpka          ) = 0._wp 
    167             z_elem_b( ji,     jpka           ) = e3t_abl( jpka ) - z_elem_a( ji,    jpka ) 
    168          END DO             
     167            z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
     168            z_elem_c( ji, jpka ) = 0._wp 
     169            z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 
     170         END DO 
    169171 
    170172         DO jtra = 1,jptq  ! loop on active tracers 
    171                 
     173 
    172174            DO jk = 3, jpkam1 
    173                DO ji = 1,jpi 
    174                   tq_abl  ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl  ( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side 
     175               !DO ji = 2, jpim1 
     176               DO ji = 1,jpi  !!GS: to be checked if needed 
     177                  tq_abl( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side 
    175178               END DO 
    176179            END DO 
    177180 
    178181            IF(jtra == jp_ta) THEN 
    179                DO ji = 1,jpi  ! boundary conditions for temperature               
    180                   zztmp1 = psen(ji, jj)  
    181                   zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )   
    182 #if defined key_si3              
    183                   zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj)  
     182               DO ji = 1,jpi  ! surface boundary condition for temperature 
     183                  zztmp1 = psen(ji, jj) 
     184                  zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 ) 
     185#if defined key_si3 
     186                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) 
    184187                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj) 
    185 #endif               
    186                   z_elem_b( ji,     2                ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1               
    187                   tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl( 2    ) * tq_abl  ( ji, jj, 2   , nt_n, jtra ) + rDt_abl * zztmp2                
     188#endif 
     189                  z_elem_b( ji,        2             ) = e3t_abl(    2 ) - z_elem_c( ji,        2             ) + rDt_abl * zztmp1 
     190                  tq_abl  ( ji, jj,    2, nt_a, jtra ) = e3t_abl(    2 ) * tq_abl  ( ji, jj,    2, nt_n, jtra ) + rDt_abl * zztmp2 
    188191                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra ) 
    189                END DO  
    190             ELSE    
    191                DO ji = 1,jpi  ! boundary conditions for humidity               
    192                   zztmp1 = pevp(ji, jj)  
    193                   zztmp2 = pevp(ji, jj) * pssq(ji, jj)    
    194 #if defined key_si3              
    195                   zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj)  
    196                   zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)     
    197 #endif    
    198                   z_elem_b( ji,     2                ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 
    199                   tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl( 2    ) * tq_abl  ( ji, jj, 2   , nt_n, jtra ) + rDt_abl * zztmp2                
     192               END DO 
     193            ELSE ! jp_qa 
     194               DO ji = 1,jpi  ! surface boundary condition for humidity 
     195                  zztmp1 = pevp(ji, jj) 
     196                  zztmp2 = pevp(ji, jj) * pssq(ji, jj) 
     197#if defined key_si3 
     198                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj) 
     199                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj) 
     200#endif 
     201                  z_elem_b( ji,     2                ) = e3t_abl(    2 ) - z_elem_c( ji,        2            ) + rDt_abl * zztmp1 
     202                  tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl(    2 ) * tq_abl  ( ji, jj,    2, nt_n, jtra ) + rDt_abl * zztmp2 
    200203                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra ) 
    201                END DO                                      
     204               END DO 
    202205            END IF 
    203206            !! 
    204207            !! Matrix inversion 
    205208            !! ---------------------------------------------------------- 
    206             DO ji = 1,jpi             
    207                zcff                       =  1._wp / z_elem_b( ji, 2 ) 
    208                zCF   (ji,   2           ) = - zcff * z_elem_c( ji, 2 ) 
    209                tq_abl(ji,jj,2,nt_a,jtra) =    zcff * tq_abl(ji,jj,2,nt_a,jtra) 
    210             END DO              
    211  
    212             DO jk = 3, jpka             
    213                DO ji = 1,jpi             
    214                   zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) ) 
     209            DO ji = 1,jpi 
     210               zcff                            =  1._wp / z_elem_b( ji, 2 ) 
     211               zCF   ( ji,     2             ) = - zcff * z_elem_c( ji, 2 ) 
     212               tq_abl( ji, jj, 2, nt_a, jtra ) =   zcff * tq_abl( ji, jj, 2, nt_a, jtra ) 
     213            END DO 
     214 
     215            DO jk = 3, jpka 
     216               DO ji = 1,jpi 
     217                  zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF( ji, jk-1 ) ) 
    215218                  zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) 
    216219                  tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk  ,nt_a,jtra)   & 
    217                   &                - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) ) 
    218                END DO 
    219             END DO 
    220 !!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...          
    221             DO jk = jpkam1,2,-1             
    222                DO ji = 1,jpi   
     220                     &           - z_elem_a(ji, jk) *  tq_abl(ji,jj,jk-1,nt_a,jtra) ) 
     221               END DO 
     222            END DO 
     223            !!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ... 
     224            DO jk = jpkam1,2,-1 
     225               DO ji = 1,jpi 
    223226                  tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) +    & 
    224227                     &                        zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra) 
    225228               END DO 
    226229            END DO 
    227           
    228          END DO   !<-- loop on tracers                     
    229          !!  
    230       !------------- 
    231       END DO             ! end outer loop  
    232       !-------------    
    233        
    234        
     230 
     231         END DO   !<-- loop on tracers 
     232         !! 
     233      !------------- 
     234      END DO             ! end outer loop 
     235      !------------- 
     236 
    235237      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    236238      !                            !  3 *** Compute Coriolis term with geostrophic guide 
    237       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   
    238       !-------------         
    239       DO jk = 2, jpka    ! outer loop 
    240       !------------- 
    241          !              
    242          ! Advance u_abl & v_abl to time n+1 
    243          DO_2D_11_11 
    244             zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl )  ! (f dt)**2 
     239      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     240      IF( SemiImp_Cor ) THEN 
     241 
     242         !------------- 
     243         DO jk = 2, jpka    ! outer loop 
     244         !------------- 
     245            ! 
     246            ! Advance u_abl & v_abl to time n+1 
     247            DO_2D_11_11 
     248               zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl )  ! (f dt)**2 
     249 
     250               u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(                                          & 
     251                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff) * u_abl( ji, jj, jk, nt_n )    & 
     252                  &                     + rDt_abl * fft_abl(ji, jj) * v_abl( ji, jj, jk, nt_n ) )  & 
     253                  &                               / (1._wp + gamma_Cor*gamma_Cor*zcff) 
    245254    
    246             u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(  & 
    247                &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n )    & 
    248                &                 +  rDt_abl * fft_abl(ji, jj) * v_abl ( ji , jj  , jk, nt_n ) )  & 
    249                &                               / (1._wp + gamma_Cor*gamma_Cor*zcff) 
    250                 
    251             v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(  & 
    252                &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n )   & 
    253                &                 -  rDt_abl * fft_abl(ji, jj) * u_abl ( ji   , jj, jk, nt_n )  ) & 
    254                &                                / (1._wp + gamma_Cor*gamma_Cor*zcff)                 
    255          END_2D 
    256          !                                    
    257       !------------- 
    258       END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 
    259       !-------------                 
    260       ! 
    261       IF( ln_geos_winds ) THEN 
    262          DO jj = 1, jpj    ! outer loop        
    263             DO jk = 1, jpka 
    264                DO ji = 1, jpi  
    265                   u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   & 
    266                      &                      - rDt_abl * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk) 
    267                   v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   & 
    268                      &                      + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk) 
    269                END DO 
    270             END DO 
    271          END DO       
    272       END IF  
    273       !-------------             
    274       ! 
    275       IF( ln_hpgls_frc ) THEN 
    276          DO jj = 1, jpj    ! outer loop        
    277             DO jk = 1, jpka 
    278                DO ji = 1, jpi  
    279                   u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk) 
    280                   v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk)   
     255               v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(                                         & 
     256                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff) * v_abl( ji, jj, jk, nt_n )    & 
     257                  &                     - rDt_abl * fft_abl(ji, jj) * u_abl( ji, jj, jk, nt_n ) )  & 
     258                  &                               / (1._wp + gamma_Cor*gamma_Cor*zcff) 
     259            END_2D 
     260            ! 
     261         !------------- 
     262         END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 
     263         !------------- 
     264         ! 
     265         IF( ln_geos_winds ) THEN 
     266            DO jj = 1, jpj    ! outer loop 
     267               DO jk = 1, jpka 
     268                  DO ji = 1, jpi 
     269                     u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   & 
     270                        &                      - rDt_abl * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk) 
     271                     v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   & 
     272                        &                      + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk) 
     273                  END DO 
     274               END DO 
     275            END DO 
     276         END IF 
     277         ! 
     278         IF( ln_hpgls_frc ) THEN 
     279            DO jj = 1, jpj    ! outer loop 
     280               DO jk = 1, jpka 
     281                  DO ji = 1, jpi 
     282                     u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk) 
     283                     v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk) 
     284                  ENDDO 
    281285               ENDDO 
    282286            ENDDO 
    283          ENDDO  
    284       END IF 
    285       
    286       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    287       !                            !  4 *** Advance u,v to time n+1  
    288       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  
    289       ! 
    290       !  Vertical diffusion for u_abl      
     287         END IF 
     288 
     289      ELSE ! SemiImp_Cor = .FALSE. 
     290 
     291         IF( ln_geos_winds ) THEN 
     292 
     293            !------------- 
     294            DO jk = 2, jpka    ! outer loop 
     295            !------------- 
     296               ! 
     297               IF( MOD( kt, 2 ) == 0 ) then 
     298                  ! Advance u_abl & v_abl to time n+1 
     299                  DO jj = 1, jpj 
     300                     DO ji = 1, jpi 
     301                        zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_n ) - pgv_dta(ji  ,jj  ,jk)  ) 
     302                        u_abl( ji, jj, jk, nt_a ) =                u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff 
     303                        zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_a ) - pgu_dta(ji  ,jj  ,jk)  ) 
     304                        v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff ) 
     305                        u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  u_abl( ji, jj, jk, nt_a ) 
     306                     END DO 
     307                  END DO 
     308               ELSE 
     309                  DO jj = 1, jpj 
     310                     DO ji = 1, jpi 
     311                        zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_n ) - pgu_dta(ji  ,jj  ,jk)  ) 
     312                        v_abl( ji, jj, jk, nt_a ) =                v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff 
     313                        zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_a ) - pgv_dta(ji  ,jj  ,jk)  ) 
     314                        u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff ) 
     315                        v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  v_abl( ji, jj, jk, nt_a ) 
     316                     END DO 
     317                  END DO 
     318               END IF 
     319               ! 
     320            !------------- 
     321            END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 
     322            !------------- 
     323 
     324         ENDIF ! ln_geos_winds 
     325 
     326      ENDIF ! SemiImp_Cor 
     327      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     328      !                            !  4 *** Advance u,v to time n+1 
     329      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     330      ! 
     331      !  Vertical diffusion for u_abl 
    291332      !------------- 
    292333      DO jj = 1, jpj    ! outer loop 
    293       !-------------     
     334      !------------- 
    294335 
    295336         DO jk = 3, jpkam1 
    296             DO ji = 1, jpi   
    297                z_elem_a( ji,     jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal 
    298                z_elem_c( ji,     jk ) = - rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal                 
    299                z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                             !       diagonal 
    300             END DO 
    301          END DO   
    302              
    303          DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi)             
     337            DO ji = 1, jpi 
     338               z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal 
     339               z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal 
     340               z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )  !       diagonal 
     341            END DO 
     342         END DO 
     343 
     344         DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi) 
    304345            !++ Surface boundary condition 
    305             z_elem_a( ji,     2    ) = 0._wp 
    306             z_elem_c( ji,     2    ) = - rDt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                        
    307             ! 
    308          zztmp1  = pcd_du(ji, jj) 
    309          zztmp2  = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) )        
    310 #if defined key_si3              
     346            z_elem_a( ji, 2 ) = 0._wp 
     347            z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 
     348            ! 
     349            zztmp1  = pcd_du(ji, jj) 
     350            zztmp2  = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) ) * rn_vfac 
     351#if defined key_si3 
    311352            zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) 
    312          zzice  = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) ) 
    313          zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 
    314 #endif            
    315          z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1          
     353            zzice  = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji, jj) ) * rn_vfac 
     354            zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 
     355#endif 
     356            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 
    316357            u_abl( ji, jj,    2, nt_a ) =      u_abl( ji, jj,    2, nt_a ) + rDt_abl * zztmp2 
    317           
    318             !++ Top Neumann B.C. 
    319             !z_elem_a( ji,     jpka ) = - 0.5_wp * rDt_abl * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka )  
    320             !z_elem_c( ji,     jpka ) = 0._wp 
    321             !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                                
    322             !++ Top Dirichlet B.C. 
    323             z_elem_a( ji,     jpka )       = 0._wp 
    324             z_elem_c( ji,     jpka )       = 0._wp 
    325             z_elem_b( ji,     jpka )       = e3t_abl( jpka ) 
    326             u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)               
    327          END DO  
     358 
     359            ! idealized test cases only 
     360            !IF( ln_topbc_neumann ) THEN 
     361            !   !++ Top Neumann B.C. 
     362            !   z_elem_a( ji,     jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
     363            !   z_elem_c( ji,     jpka ) = 0._wp 
     364            !   z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka ) 
     365            !   !u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * u_abl   ( ji, jj, jpka, nt_a ) 
     366            !ELSE 
     367               !++ Top Dirichlet B.C. 
     368               z_elem_a( ji,     jpka )       = 0._wp 
     369               z_elem_c( ji,     jpka )       = 0._wp 
     370               z_elem_b( ji,     jpka )       = e3t_abl( jpka ) 
     371               u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk) 
     372            !ENDIF 
     373 
     374         END DO 
    328375         !! 
    329376         !! Matrix inversion 
    330377         !! ---------------------------------------------------------- 
    331          DO ji = 2, jpi           
     378         !DO ji = 2, jpi 
     379         DO ji = 1, jpi  !!GS: TBI 
    332380            zcff                 =   1._wp / z_elem_b( ji, 2 ) 
    333             zCF   (ji,   2     ) =  - zcff * z_elem_c( ji, 2 )  
     381            zCF   (ji,   2     ) =  - zcff * z_elem_c( ji, 2 ) 
    334382            u_abl (ji,jj,2,nt_a) =    zcff * u_abl(ji,jj,2,nt_a) 
    335          END DO              
    336  
    337          DO jk = 3, jpka             
    338             DO ji = 2, jpi               
     383         END DO 
     384 
     385         DO jk = 3, jpka 
     386            DO ji = 2, jpi 
    339387               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) ) 
    340388               zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) 
     
    343391            END DO 
    344392         END DO 
    345              
    346          DO jk = jpkam1,2,-1             
     393 
     394         DO jk = jpkam1,2,-1 
    347395            DO ji = 2, jpi 
    348396               u_abl(ji,jj,jk,nt_a) = u_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * u_abl(ji,jj,jk+1,nt_a) 
    349397            END DO 
    350398         END DO 
    351                                            
    352       !------------- 
    353       END DO             ! end outer loop  
    354       !-------------    
    355  
    356       ! 
    357       !  Vertical diffusion for v_abl      
     399 
     400      !------------- 
     401      END DO             ! end outer loop 
     402      !------------- 
     403 
     404      ! 
     405      !  Vertical diffusion for v_abl 
    358406      !------------- 
    359407      DO jj = 2, jpj   ! outer loop 
    360       !-------------     
     408      !------------- 
    361409         ! 
    362410         DO jk = 3, jpkam1 
    363             DO ji = 1, jpi    
    364                z_elem_a( ji,     jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal 
    365                z_elem_c( ji,     jk ) = -rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal               
    366                z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal 
    367             END DO 
    368          END DO 
    369  
    370          DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)           
     411            DO ji = 1, jpi 
     412               z_elem_a( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal 
     413               z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal 
     414               z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal 
     415            END DO 
     416         END DO 
     417 
     418         DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj) 
    371419            !++ Surface boundary condition 
    372             z_elem_a( ji,     2    ) = 0._wp 
    373             z_elem_c( ji,     2    ) = - rDt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )         
    374             ! 
    375          zztmp1 = pcd_du(ji, jj) 
    376          zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) )   
    377 #if defined key_si3              
     420            z_elem_a( ji, 2 ) = 0._wp 
     421            z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 
     422            ! 
     423            zztmp1 = pcd_du(ji, jj) 
     424            zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) ) * rn_vfac 
     425#if defined key_si3 
    378426            zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) 
    379          zzice  = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji,jj-1) ) 
    380          zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 
    381 #endif          
    382             z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1  
     427            zzice  = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji, jj-1) ) * rn_vfac 
     428            zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 
     429#endif 
     430            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 
    383431            v_abl( ji, jj,    2, nt_a ) =         v_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2 
    384             !++ Top Neumann B.C.             
    385             !z_elem_a( ji,     jpka ) = -rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )  
    386             !z_elem_c( ji,     jpka ) = 0._wp 
    387             !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                                
    388             !++ Top Dirichlet B.C. 
    389             z_elem_a( ji,     jpka )       = 0._wp 
    390             z_elem_c( ji,     jpka )       = 0._wp 
    391             z_elem_b( ji,     jpka )       = e3t_abl( jpka ) 
    392             v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)              
    393          END DO  
     432 
     433            ! idealized test cases only 
     434            !IF( ln_topbc_neumann ) THEN 
     435            !   !++ Top Neumann B.C. 
     436            !   z_elem_a( ji,     jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
     437            !   z_elem_c( ji,     jpka ) = 0._wp 
     438            !   z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka ) 
     439            !   !v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * v_abl   ( ji, jj, jpka, nt_a ) 
     440            !ELSE 
     441               !++ Top Dirichlet B.C. 
     442               z_elem_a( ji,     jpka )       = 0._wp 
     443               z_elem_c( ji,     jpka )       = 0._wp 
     444               z_elem_b( ji,     jpka )       = e3t_abl( jpka ) 
     445               v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk) 
     446            !ENDIF 
     447 
     448         END DO 
    394449         !! 
    395450         !! Matrix inversion 
    396451         !! ---------------------------------------------------------- 
    397          DO ji = 1, jpi               
     452         DO ji = 1, jpi 
    398453            zcff                 =  1._wp / z_elem_b( ji, 2 ) 
    399             zCF   (ji,   2     ) =   - zcff * z_elem_c( ji,     2       )  
    400             v_abl (ji,jj,2,nt_a) =     zcff * v_abl   ( ji, jj, 2, nt_a )  
    401          END DO              
    402  
    403          DO jk = 3, jpka             
    404             DO ji = 1, jpi                    
     454            zCF   (ji,   2     ) =   - zcff * z_elem_c( ji,     2       ) 
     455            v_abl (ji,jj,2,nt_a) =     zcff * v_abl   ( ji, jj, 2, nt_a ) 
     456         END DO 
     457 
     458         DO jk = 3, jpka 
     459            DO ji = 1, jpi 
    405460               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) ) 
    406461               zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) 
     
    409464            END DO 
    410465         END DO 
    411              
    412          DO jk = jpkam1,2,-1             
    413             DO ji = 1, jpi     
     466 
     467         DO jk = jpkam1,2,-1 
     468            DO ji = 1, jpi 
    414469               v_abl(ji,jj,jk,nt_a) = v_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * v_abl(ji,jj,jk+1,nt_a) 
    415470            END DO 
    416471         END DO 
    417          !                                           
    418       !------------- 
    419       END DO             ! end outer loop  
    420       !-------------    
    421      
    422       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    423       !                            !  5 *** Apply nudging on the dynamics and the tracers  
    424       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  
    425       z_cft(:,:,:) = 0._wp       
    426        
     472         ! 
     473      !------------- 
     474      END DO             ! end outer loop 
     475      !------------- 
     476 
     477      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     478      !                            !  5 *** Apply nudging on the dynamics and the tracers 
     479      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     480 
    427481      IF( nn_dyn_restore > 0  ) THEN 
    428          !-------------  
     482         !------------- 
    429483         DO jk = 2, jpka    ! outer loop 
    430          !-------------        
     484         !------------- 
    431485            DO_2D_01_01 
    432486               zcff1 = pblh( ji, jj ) 
    433                zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )                         
    434                zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )  
     487               zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) ) 
     488               zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
    435489               zmsk  = msk_abl(ji,jj) 
    436490               zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2   & 
    437491                  &  + jp_alp1_dyn * zsig    + jp_alp0_dyn 
    438492               zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points 
    439                                                              ! rn_Dt = rDt_abl / nn_fsbc                           
     493                                                             ! rn_Dt = rDt_abl / nn_fsbc 
    440494               zcff  = zcff * rest_eq(ji,jj) 
    441                z_cft( ji, jj, jk ) = zcff 
    442495               u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  u_abl( ji, jj, jk, nt_a )   & 
    443                   &                               + zcff   * pu_dta( ji, jj, jk       )                       
     496                  &                               + zcff   * pu_dta( ji, jj, jk       ) 
    444497               v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  v_abl( ji, jj, jk, nt_a )   & 
    445498                  &                               + zcff   * pv_dta( ji, jj, jk       ) 
     
    447500         !------------- 
    448501         END DO             ! end outer loop 
    449          !-------------                
     502         !------------- 
    450503      END IF 
    451504 
    452       !-------------  
     505      !------------- 
    453506      DO jk = 2, jpka    ! outer loop 
    454       !-------------        
     507      !------------- 
    455508         DO_2D_11_11 
    456509            zcff1 = pblh( ji, jj ) 
    457510            zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) ) 
    458             zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )  
     511            zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
    459512            zmsk  = msk_abl(ji,jj) 
    460513            zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2   & 
    461514               &  + jp_alp1_tra * zsig    + jp_alp0_tra 
    462515            zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points 
    463                                                           ! rn_Dt = rDt_abl / nn_fsbc                           
    464             !z_cft( ji, jj, jk ) = zcff 
     516                                                          ! rn_Dt = rDt_abl / nn_fsbc 
    465517            tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta )   & 
    466518               &                                       + zcff   * pt_dta( ji, jj, jk ) 
    467              
     519 
    468520            tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa )   & 
    469521               &                                       + zcff   * pq_dta( ji, jj, jk ) 
    470              
     522 
    471523         END_2D 
    472524      !------------- 
    473525      END DO             ! end outer loop 
    474       !-------------       
    475       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    476       !                            !  6 *** MPI exchanges   
    477       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    478       ! 
    479       CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a      ), 'T', -1. ) 
     526      !------------- 
     527      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     528      !                            !  6 *** MPI exchanges 
     529      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     530      ! 
     531      CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a)      , 'T', -1.                            ) 
    480532      CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T',  1., tq_abl(:,:,:,nt_a,jp_qa), 'T',  1., kfillmode = jpfillnothing )   ! ++++ this should not be needed... 
    481533      ! 
    482       ! first ABL level 
     534#if defined key_iomput 
     535      ! 2D & first ABL level 
     536      IF ( iom_use("pblh"   ) ) CALL iom_put (    "pblh",    pblh(:,:             ) ) 
    483537      IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl",   u_abl(:,:,2,nt_a      ) ) 
    484538      IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl",   v_abl(:,:,2,nt_a      ) ) 
     
    489543      IF ( iom_use("tz1_dta") ) CALL iom_put ( "tz1_dta",  pt_dta(:,:,2           ) ) 
    490544      IF ( iom_use("qz1_dta") ) CALL iom_put ( "qz1_dta",  pq_dta(:,:,2           ) ) 
    491       ! all ABL levels 
    492       IF ( iom_use("u_abl"  ) ) CALL iom_put ( "u_abl"  ,   u_abl(:,:,2:jpka,nt_a      ) ) 
    493       IF ( iom_use("v_abl"  ) ) CALL iom_put ( "v_abl"  ,   v_abl(:,:,2:jpka,nt_a      ) ) 
    494       IF ( iom_use("t_abl"  ) ) CALL iom_put ( "t_abl"  ,  tq_abl(:,:,2:jpka,nt_a,jp_ta) ) 
    495       IF ( iom_use("q_abl"  ) ) CALL iom_put ( "q_abl"  ,  tq_abl(:,:,2:jpka,nt_a,jp_qa) ) 
    496       IF ( iom_use("tke_abl") ) CALL iom_put ( "tke_abl", tke_abl(:,:,2:jpka,nt_a      ) ) 
    497       IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka           ) ) 
    498       IF ( iom_use("avt_abl") ) CALL iom_put ( "avt_abl", avm_abl(:,:,2:jpka           ) ) 
    499       IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxl_abl(:,:,2:jpka           ) ) 
    500       IF ( iom_use("pblh"   ) ) CALL iom_put (  "pblh"  ,    pblh(:,:                  ) ) 
    501       ! debug (to be removed) 
     545      ! debug 2D 
     546      IF( ln_geos_winds ) THEN 
     547         IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2) ) 
     548         IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2) ) 
     549      END IF 
     550      IF( ln_hpgls_frc ) THEN 
     551         IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo",  pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) ) 
     552         IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) ) 
     553      END IF 
     554      ! 3D (all ABL levels) 
     555      IF ( iom_use("u_abl"   ) ) CALL iom_put ( "u_abl"   ,    u_abl(:,:,2:jpka,nt_a      ) ) 
     556      IF ( iom_use("v_abl"   ) ) CALL iom_put ( "v_abl"   ,    v_abl(:,:,2:jpka,nt_a      ) ) 
     557      IF ( iom_use("t_abl"   ) ) CALL iom_put ( "t_abl"   ,   tq_abl(:,:,2:jpka,nt_a,jp_ta) ) 
     558      IF ( iom_use("q_abl"   ) ) CALL iom_put ( "q_abl"   ,   tq_abl(:,:,2:jpka,nt_a,jp_qa) ) 
     559      IF ( iom_use("tke_abl" ) ) CALL iom_put ( "tke_abl" ,  tke_abl(:,:,2:jpka,nt_a      ) ) 
     560      IF ( iom_use("avm_abl" ) ) CALL iom_put ( "avm_abl" ,  avm_abl(:,:,2:jpka           ) ) 
     561      IF ( iom_use("avt_abl" ) ) CALL iom_put ( "avt_abl" ,  avt_abl(:,:,2:jpka           ) ) 
     562      IF ( iom_use("mxlm_abl") ) CALL iom_put ( "mxlm_abl", mxlm_abl(:,:,2:jpka           ) ) 
     563      IF ( iom_use("mxld_abl") ) CALL iom_put ( "mxld_abl", mxld_abl(:,:,2:jpka           ) ) 
     564      ! debug 3D 
    502565      IF ( iom_use("u_dta") ) CALL iom_put ( "u_dta",  pu_dta(:,:,2:jpka) ) 
    503566      IF ( iom_use("v_dta") ) CALL iom_put ( "v_dta",  pv_dta(:,:,2:jpka) ) 
    504567      IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta",  pt_dta(:,:,2:jpka) ) 
    505568      IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta",  pq_dta(:,:,2:jpka) ) 
    506       IF ( iom_use("coeft") ) CALL iom_put ( "coeft",   z_cft(:,:,2:jpka) ) 
    507569      IF( ln_geos_winds ) THEN 
    508          IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2           ) ) 
    509          IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2           ) ) 
     570         IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", pgu_dta(:,:,2:jpka) ) 
     571         IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", pgv_dta(:,:,2:jpka) ) 
    510572      END IF 
    511573      IF( ln_hpgls_frc ) THEN 
    512          IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo",  pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp)  ) 
    513          IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp)  ) 
     574         IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo",  pgu_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) ) 
     575         IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", -pgv_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) ) 
    514576      END IF 
    515       ! 
     577#endif 
    516578      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    517579      !                            !  7 *** Finalize flux computation 
    518       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  
    519  
     580      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     581      ! 
    520582      DO_2D_11_11 
    521          ztemp             = tq_abl  ( ji, jj, 2, nt_a, jp_ta )  
    522          zhumi             = tq_abl  ( ji, jj, 2, nt_a, jp_qa )  
    523          !zcff              = pslp_dta( ji, jj ) /   &              !<-- At this point ztemp and zhumi should not be zero ... 
    524          !   &                        (  R_dry*ztemp * ( 1._wp + rctv0*zhumi )  ) 
    525          zcff              = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 
    526          psen ( ji, jj )   =      cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) 
    527          pevp ( ji, jj )   = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) ) 
    528          rhoa( ji, jj )   = zcff               
     583         ztemp          =  tq_abl( ji, jj, 2, nt_a, jp_ta ) 
     584         zhumi          =  tq_abl( ji, jj, 2, nt_a, jp_qa ) 
     585         zcff           = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 
     586         psen( ji, jj ) =    - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )   !GS: negative sign to respect aerobulk convention 
     587         pevp( ji, jj ) = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) ) 
     588         rhoa( ji, jj ) = zcff 
    529589      END_2D 
    530        
     590 
    531591      DO_2D_01_01 
     592<<<<<<< .working 
    532593         zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * ( pssu(ji  ,jj) + pssu(ji-1,jj) )   
    533594         zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * ( pssv(ji,jj  ) + pssv(ji,jj-1) )  
     595||||||| .merge-left.r13136 
     596         zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji  ,jj) + pssu(ji-1,jj) )   
     597         zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj  ) + pssv(ji,jj-1) )  
     598======= 
     599         zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji,jj) + pssu(ji-1,jj  ) ) 
     600         zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj) + pssv(ji  ,jj-1) ) 
     601>>>>>>> .merge-right.r13211 
    534602      END_2D 
    535       !  
     603      ! 
    536604      CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. ) 
    537605      ! 
     
    539607      DO_2D_11_11 
    540608         zcff          = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   & 
    541             &                 + zwnd_j(ji,jj) * zwnd_j(ji,jj)  )  ! * msk_abl(ji,jj) 
     609            &                 + zwnd_j(ji,jj) * zwnd_j(ji,jj) )   ! * msk_abl(ji,jj) 
    542610         zztmp         = rhoa(ji,jj) * pcd_du(ji,jj) 
    543           
     611 
    544612         pwndm (ji,jj) =         zcff 
    545613         ptaum (ji,jj) = zztmp * zcff 
     
    564632 
    565633      IF(sn_cfctl%l_prtctl) THEN 
    566          CALL prt_ctl( tab2d_1=pwndm  , clinfo1=' abl_stp: wndm   : ' ) 
    567          CALL prt_ctl( tab2d_1=ptaui  , clinfo1=' abl_stp: utau   : ' ) 
    568          CALL prt_ctl( tab2d_2=ptauj  , clinfo2=          'vtau   : ' ) 
     634         CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau   : ', mask1=umask,   & 
     635            &          tab2d_2=ptauj , clinfo2='          vtau   : ', mask2=vmask ) 
     636         CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm   : ' ) 
    569637      ENDIF 
    570638 
    571639#if defined key_si3 
     640<<<<<<< .working 
    572641      ! ------------------------------------------------------------ ! 
    573642      !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     
    583652      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   & 
    584653         &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' ) 
     654||||||| .merge-left.r13136 
     655         ! ------------------------------------------------------------ ! 
     656         !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     657         ! ------------------------------------------------------------ ! 
     658         DO_2D_00_00 
     659             
     660            zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 
     661            zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 
     662    
     663            ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             & 
     664               &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          & 
     665               &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 
     666            ptauj_ice(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             & 
     667               &                      +    rhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          & 
     668               &         * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 
     669         END_2D 
     670         CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. ) 
     671         ! 
     672         IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   & 
     673            &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' ) 
     674======= 
     675      ! ------------------------------------------------------------ ! 
     676      !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     677      ! ------------------------------------------------------------ ! 
     678      DO_2D_00_00 
     679 
     680         zztmp1 = 0.5_wp * ( u_abl(ji+1,jj  ,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 
     681         zztmp2 = 0.5_wp * ( v_abl(ji  ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 
     682 
     683         ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             & 
     684            &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          & 
     685            &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 
     686         ptauj_ice(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             & 
     687            &                      +    rhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          & 
     688            &         * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 
     689      END_2D 
     690      CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. ) 
     691      ! 
     692      IF(sn_cfctl%l_prtctl) THEN 
     693         CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=umask,   & 
     694            &          tab2d_2=ptauj_ice , clinfo2='          vtau_ice : ', mask2=vmask ) 
     695      END IF 
     696>>>>>>> .merge-right.r13211 
    585697#endif 
    586698      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    615727      !!                (= Kz dz[Ub] * dz[Un] ) 
    616728      !! --------------------------------------------------------------------- 
    617       INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn  
     729      INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn 
    618730      INTEGER, DIMENSION(1:jpi          )     ::   ikbl 
    619731      REAL(wp)                                ::   zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv 
    620       REAL(wp)                                ::   zdU, zdV, zcff1,zshear,zbuoy,zsig, zustar2 
    621       REAL(wp)                                ::   zdU2,zdV2       
    622       REAL(wp)                                ::   zwndi,zwndj 
     732      REAL(wp)                                ::   zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2 
     733      REAL(wp)                                ::   zdU2, zdV2, zbuoy1, zbuoy2    ! zbuoy for BL89 
     734      REAL(wp)                                ::   zwndi, zwndj 
    623735      REAL(wp), DIMENSION(1:jpi,      1:jpka) ::   zsh2 
    624736      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) ::   zbn2 
    625       REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF        
     737      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF 
    626738      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_a 
    627739      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_b 
     
    629741      LOGICAL                                 ::   ln_Patankar    = .FALSE. 
    630742      LOGICAL                                 ::   ln_dumpvar     = .FALSE. 
    631       LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl  
     743      LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl 
    632744      ! 
    633745      tind  = nt_n 
     
    641753      !------------- 
    642754         ! 
    643          ! Compute vertical shear          
     755         ! Compute vertical shear 
    644756         DO jk = 2, jpkam1 
    645             DO ji = 1,jpi   
    646                zcff        = 1.0_wp / e3w_abl( jk )**2                 
    647                zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2                            
     757            DO ji = 1, jpi 
     758               zcff        = 1.0_wp / e3w_abl( jk )**2 
     759               zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2 
    648760               zdV         = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2 
    649                zsh2(ji,jk) = zdU+zdV 
    650             END DO 
    651          END DO    
     761               zsh2(ji,jk) = zdU+zdV   !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 ) 
     762            END DO 
     763         END DO 
    652764         ! 
    653765         ! Compute brunt-vaisala frequency 
    654766         DO jk = 2, jpkam1 
    655             DO ji = 1,jpi  
    656                zcff  = grav * itvref / e3w_abl( jk )   
     767            DO ji = 1,jpi 
     768               zcff  = grav * itvref / e3w_abl( jk ) 
    657769               zcff1 =  tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk  , tind, jp_ta) 
    658770               zcff2 =  tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa)        & 
     
    660772               zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 )  !<-- zbn2 defined on (2,jpi) 
    661773            END DO 
    662          END DO  
     774         END DO 
    663775         ! 
    664776         ! Terms for the tridiagonal problem 
    665777         DO jk = 2, jpkam1 
    666             DO ji = 1,jpi  
    667                zshear       =                 zsh2( ji,     jk )   ! zsh2 is already multiplied by Avm_abl at this point 
    668                zsh2(ji,jk)  = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation           
    669                zbuoy        = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )  
    670                 
    671                z_elem_a( ji,     jk )   = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk   ) ! lower-diagonal 
    672                z_elem_c( ji,     jk )   = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal           
     778            DO ji = 1, jpi 
     779               zshear      = zsh2( ji, jk )                           ! zsh2 is already multiplied by Avm_abl at this point 
     780               zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation 
     781               zbuoy       = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 
     782 
     783               z_elem_a( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk ) + Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk   ) ! lower-diagonal 
     784               z_elem_c( ji, jk ) = - 0.5_wp * rDt_abl * rn_Sch * ( Avm_abl( ji, jj, jk ) + Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal 
    673785               IF( (zbuoy + zshear) .gt. 0.) THEN    ! Patankar trick to avoid negative values of TKE 
    674                   z_elem_b( ji,     jk )  = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
    675                      &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)     ! diagonal        
    676                   tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )             ! right-hand-side 
     786                  z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
     787                     &               + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk)    ! diagonal 
     788                  tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )       ! right-hand-side 
    677789               ELSE 
    678                   z_elem_b( ji,     jk )  = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
    679                      &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)   &  ! diagonal     
    680                      &                     - e3w_abl(jk) * rDt_abl * zbuoy    
    681                   tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl *  zshear )             ! right-hand-side                      
     790                  z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
     791                     &               + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk)   &  ! diagonal 
     792                     &               - e3w_abl(jk) * rDt_abl * zbuoy 
     793                  tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl *  zshear )                    ! right-hand-side 
    682794               END IF 
    683795            END DO 
    684          END DO   
    685                              
    686          DO ji = 1,jpi    ! vector opt.             
    687             zesrf   =  MAX( 4.63_wp * ustar2(ji,jj), tke_min )            
    688             zetop   = tke_min              
    689             z_elem_a ( ji,     1    ) = 0._wp; z_elem_c ( ji,     1    ) = 0._wp; z_elem_b ( ji,     1    ) = 1._wp                                                    
    690             z_elem_a ( ji,     jpka ) = 0._wp; z_elem_c ( ji,     jpka ) = 0._wp; z_elem_b ( ji,     jpka ) = 1._wp             
    691             tke_abl( ji, jj,    1, nt_a ) = zesrf 
    692             tke_abl( ji, jj, jpka, nt_a ) = zetop  
    693             zbn2(ji,jj,   1) = zbn2( ji,jj,   2)  
    694             zsh2(ji,      1) = zsh2( ji,      2)                                       
    695             zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1)  
    696             zsh2(ji,   jpka) = zsh2( ji  , jpkam1) 
    697          END DO         
     796         END DO 
     797 
     798         DO ji = 1,jpi    ! vector opt. 
     799            zesrf = MAX( rn_Esfc * ustar2(ji,jj), tke_min ) 
     800            zetop = tke_min 
     801 
     802            z_elem_a ( ji,     1       ) = 0._wp 
     803            z_elem_c ( ji,     1       ) = 0._wp 
     804            z_elem_b ( ji,     1       ) = 1._wp 
     805            tke_abl  ( ji, jj, 1, nt_a ) = zesrf 
     806 
     807            !++ Top Neumann B.C. 
     808            !z_elem_a ( ji,     jpka       ) = - 0.5 * rDt_abl * rn_Sch * (Avm_abl(ji,jj, jpka-1 )+Avm_abl(ji,jj, jpka ))  / e3t_abl( jpka ) 
     809            !z_elem_c ( ji,     jpka       ) = 0._wp 
     810            !z_elem_b ( ji,     jpka       ) = e3w_abl(jpka) - z_elem_a(ji, jpka ) 
     811            !tke_abl  ( ji, jj, jpka, nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n ) 
     812 
     813            !++ Top Dirichlet B.C. 
     814            z_elem_a ( ji,     jpka       ) = 0._wp 
     815            z_elem_c ( ji,     jpka       ) = 0._wp 
     816            z_elem_b ( ji,     jpka       ) = 1._wp 
     817            tke_abl  ( ji, jj, jpka, nt_a ) = zetop 
     818 
     819            zbn2 ( ji, jj,    1 ) = zbn2 ( ji, jj,      2 ) 
     820            zsh2 ( ji,        1 ) = zsh2 ( ji,          2 ) 
     821            zbn2 ( ji, jj, jpka ) = zbn2 ( ji, jj, jpkam1 ) 
     822            zsh2 ( ji,     jpka ) = zsh2 ( ji    , jpkam1 ) 
     823         END DO 
    698824         !! 
    699825         !! Matrix inversion 
     
    701827         DO ji = 1,jpi 
    702828            zcff                  =  1._wp / z_elem_b( ji, 1 ) 
    703             zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       )  
    704             tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a )  
    705          END DO              
    706  
    707          DO jk = 2, jpka             
     829            zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       ) 
     830            tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a ) 
     831         END DO 
     832 
     833         DO jk = 2, jpka 
    708834            DO ji = 1,jpi 
    709835               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) ) 
     
    713839            END DO 
    714840         END DO 
    715              
    716          DO jk = jpkam1,1,-1             
     841 
     842         DO jk = jpkam1,1,-1 
    717843            DO ji = 1,jpi 
    718844               tke_abl(ji,jj,jk,nt_a) = tke_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * tke_abl(ji,jj,jk+1,nt_a) 
    719845            END DO 
    720846         END DO 
    721           
    722 !!FL should not be needed because of Patankar procedure   
     847 
     848!!FL should not be needed because of Patankar procedure 
    723849         tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min ) 
    724850 
     
    726852         !! Diagnose PBL height 
    727853         !! ---------------------------------------------------------- 
    728           
    729           
    730          !                                                        
     854 
     855 
     856         ! 
    731857         ! arrays zRH, zFC and zCF are available at this point 
    732858         ! and zFC(:, 1 ) = 0. 
    733859         ! diagnose PBL height based on zsh2 and zbn2 
    734860         zFC (  :  ,1) = 0._wp 
    735          ikbl( 1:jpi ) = 0  
    736           
     861         ikbl( 1:jpi ) = 0 
     862 
    737863         DO jk = 2,jpka 
    738             DO ji = 1, jpi  
     864            DO ji = 1, jpi 
    739865               zcff  = ghw_abl( jk-1 ) 
    740866               zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) ) 
     
    762888            ELSE 
    763889               pblh( ji, jj ) = ghw_abl(jpka) 
    764             END IF  
    765          END DO 
    766       !-------------       
    767       END DO       
    768       !------------- 
    769       ! 
    770       ! Optional : could add pblh smoothing if pblh is noisy horizontally ...  
     890            END IF 
     891         END DO 
     892      !------------- 
     893      END DO 
     894      !------------- 
     895      ! 
     896      ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 
    771897      IF(ln_smth_pblh) THEN 
    772          CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) 
     898         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) !, kfillmode = jpfillnothing) 
    773899         CALL smooth_pblh( pblh, msk_abl ) 
    774          CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)    
     900         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) !, kfillmode = jpfillnothing) 
    775901      ENDIF 
    776902      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    780906      SELECT CASE ( nn_amxl ) 
    781907      ! 
    782       CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom          
    783 #   define zlup zRH       
    784 #   define zldw zFC  
     908      CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom 
     909#   define zlup zRH 
     910#   define zldw zFC 
    785911         DO jj = 1, jpj     ! outer loop 
    786912            ! 
    787913            DO ji = 1, jpi 
    788                mxl_abl ( ji, jj,    1 )  = 0._wp 
    789                mxl_abl ( ji, jj, jpka )  = mxl_min 
    790                zldw( ji,        1 )  = 0._wp               
    791                zlup( ji,     jpka )  = 0._wp 
    792             END DO 
    793             ! 
    794             DO jk = 2, jpkam1    
    795                DO ji = 1, jpi   
    796                   zbuoy             = MAX( zbn2(ji, jj, jk), rsmall )   
    797                   mxl_abl( ji, jj, jk ) = MAX( mxl_min,  & 
    798                      &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy )   )  
    799                END DO 
    800             END DO    
     914               mxld_abl( ji, jj,    1 ) = mxl_min 
     915               mxld_abl( ji, jj, jpka ) = mxl_min 
     916               mxlm_abl( ji, jj,    1 ) = mxl_min 
     917               mxlm_abl( ji, jj, jpka ) = mxl_min 
     918               zldw    ( ji,        1 ) = zrough(ji,jj) * rn_Lsfc 
     919               zlup    ( ji,     jpka ) = mxl_min 
     920            END DO 
     921            ! 
     922            DO jk = 2, jpkam1 
     923               DO ji = 1, jpi 
     924                  zbuoy = MAX( zbn2(ji, jj, jk), rsmall ) 
     925                  mxlm_abl( ji, jj, jk ) = MAX( mxl_min,  & 
     926                     &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) ) 
     927               END DO 
     928            END DO 
    801929            ! 
    802930            ! Limit mxl 
    803             DO jk = jpkam1,1,-1    
    804                DO ji = 1, jpi  
    805                   zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl_abl(ji, jj, jk) ) 
    806                END DO 
    807             END DO    
     931            DO jk = jpkam1,1,-1 
     932               DO ji = 1, jpi 
     933                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 
     934               END DO 
     935            END DO 
    808936            ! 
    809937            DO jk = 2, jpka 
    810                DO ji = 1, jpi   
    811                   zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) )    
    812                END DO 
    813             END DO  
     938               DO ji = 1, jpi 
     939                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 
     940               END DO 
     941            END DO 
     942            ! 
     943!            DO jk = 1, jpka 
     944!               DO ji = 1, jpi 
     945!                  mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     946!                  mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ),  zlup( ji, jk ) ) 
     947!               END DO 
     948!            END DO 
    814949            ! 
    815950            DO jk = 1, jpka 
    816951               DO ji = 1, jpi 
    817                   mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )    
    818                END DO 
    819             END DO                          
    820             ! 
    821          END DO 
    822 #   undef zlup       
    823 #   undef zldw   
    824          ! 
    825          ! 
    826       CASE ( 1 )             ! length-scale computed as the distance to the PBL height 
    827          DO jj = 1,jpj      ! outer loop 
    828             ! 
    829             DO ji = 1, jpi   ! vector opt. 
    830                zcff = 1._wp / pblh( ji, jj )     ! inverse of hbl               
    831                DO jk = 1, jpka               
    832                   zsig  = MIN( zcff * ghw_abl( jk ), 1. )     
    833                   zcff1 = pblh( ji, jj )                  
    834                   mxl_abl( ji, jj, jk ) =  mxl_min                           & 
    835                      &  +  zsig         * (  amx1*zcff1 + bmx1*mxl_min ) & 
    836                      &  +  zsig *  zsig * (  amx2*zcff1 + bmx2*mxl_min ) & 
    837                      &  +  zsig**3      * (  amx3*zcff1 + bmx3*mxl_min ) & 
    838                      &  +  zsig**4      * (  amx4*zcff1 + bmx4*mxl_min ) & 
    839                      &  +  zsig**5      * (  amx5*zcff1 + bmx5*mxl_min ) 
    840                END DO 
    841             END DO  
    842             !             
    843          END DO             
     952!                  zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp) 
     953                  zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     954                  mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 
     955                  mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min ) 
     956               END DO 
     957            END DO 
     958            ! 
     959         END DO 
     960#   undef zlup 
     961#   undef zldw 
     962         ! 
     963         ! 
     964      CASE ( 1 )           ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom 
     965#   define zlup zRH 
     966#   define zldw zFC 
     967         DO jj = 1, jpj     ! outer loop 
     968            ! 
     969            DO jk = 2, jpkam1 
     970               DO ji = 1,jpi 
     971                              zcff        = 1.0_wp / e3w_abl( jk )**2 
     972                  zdU         = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2 
     973                  zdV         = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2 
     974                  zsh2(ji,jk) = SQRT(zdU+zdV)   !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 
     975                           ENDDO 
     976                        ENDDO 
     977                        ! 
     978            DO ji = 1, jpi 
     979               zcff                      = zrough(ji,jj) * rn_Lsfc 
     980               mxld_abl ( ji, jj,    1 ) = zcff 
     981               mxld_abl ( ji, jj, jpka ) = mxl_min 
     982               mxlm_abl ( ji, jj,    1 ) = zcff 
     983               mxlm_abl ( ji, jj, jpka ) = mxl_min 
     984               zldw     ( ji,        1 ) = zcff 
     985               zlup     ( ji,     jpka ) = mxl_min 
     986            END DO 
     987            ! 
     988            DO jk = 2, jpkam1 
     989               DO ji = 1, jpi 
     990                  zbuoy    = MAX( zbn2(ji, jj, jk), rsmall ) 
     991                  zcff     = 2.*SQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rod*zsh2(ji,jk) & 
     992                                &             + SQRT( rn_Rod*rn_Rod*zsh2(ji,jk)*zsh2(ji,jk)+2.*zbuoy ) ) 
     993                                  mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff ) 
     994               END DO 
     995            END DO 
     996            ! 
     997            ! Limit mxl 
     998            DO jk = jpkam1,1,-1 
     999               DO ji = 1, jpi 
     1000                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 
     1001               END DO 
     1002            END DO 
     1003            ! 
     1004            DO jk = 2, jpka 
     1005               DO ji = 1, jpi 
     1006                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 
     1007               END DO 
     1008            END DO 
     1009            ! 
     1010            DO jk = 1, jpka 
     1011               DO ji = 1, jpi 
     1012                  !mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     1013                  !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp) 
     1014                  zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     1015                  mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 
     1016                  !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) ) 
     1017                  mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min ) 
     1018               END DO 
     1019            END DO 
     1020 ! 
     1021         END DO 
     1022#   undef zlup 
     1023#   undef zldw 
    8441024         ! 
    8451025      CASE ( 2 )           ! Bougeault & Lacarrere 89 length-scale 
    8461026         ! 
    847 #   define zlup zRH       
    848 #   define zldw zFC            
     1027#   define zlup zRH 
     1028#   define zldw zFC 
    8491029! zCF is used for matrix inversion 
    850 !              
     1030! 
    8511031       DO jj = 1, jpj      ! outer loop 
    852           
    853          DO ji = 1, jpi            
    854             zlup( ji,    1 ) = mxl_min 
    855             zldw( ji,    1 ) = mxl_min                
     1032 
     1033         DO ji = 1, jpi 
     1034            zcff             = zrough(ji,jj) * rn_Lsfc 
     1035            zlup( ji,    1 ) = zcff 
     1036            zldw( ji,    1 ) = zcff 
    8561037            zlup( ji, jpka ) = mxl_min 
    857             zldw( ji, jpka ) = mxl_min                 
    858          END DO             
    859           
     1038            zldw( ji, jpka ) = mxl_min 
     1039         END DO 
     1040 
    8601041         DO jk = 2,jpka-1 
    8611042            DO ji = 1, jpi 
    8621043               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk) 
    863                zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)         
    864             END DO 
    865          END DO          
     1044               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1) 
     1045            END DO 
     1046         END DO 
    8661047         !! 
    8671048         !! BL89 search for lup 
    868          !! ----------------------------------------------------------            
    869          DO jk=2,jpka-1  
     1049         !! ---------------------------------------------------------- 
     1050         DO jk=2,jpka-1 
    8701051            ! 
    8711052            DO ji = 1, jpi 
     
    8731054               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a ) 
    8741055               ln_foundl(ji ) = .false. 
    875             END DO    
    876             !            
     1056            END DO 
     1057            ! 
    8771058            DO jkup=jk+1,jpka-1 
    8781059               DO ji = 1, jpi 
     1060                  zbuoy1 = MAX( zbn2(ji,jj,jkup  ), rsmall ) 
     1061                  zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall ) 
    8791062                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * & 
    880                      &               ( zbn2(ji,jj,jkup  )*(ghw_abl(jkup  )-ghw_abl(jk)) & 
    881                      &               + zbn2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) ) 
     1063                     &               ( zbuoy1*(ghw_abl(jkup  )-ghw_abl(jk)) & 
     1064                     &               + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) 
    8821065                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 
    8831066                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk) 
    884                      zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)                    
     1067                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 
    8851068                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   & 
    886                         &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) )  
    887                      zlup(ji,jk) = zcff                 
     1069                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) ) 
     1070                     zlup(ji,jk) = zcff 
     1071                     zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk) 
    8881072                     ln_foundl(ji) = .true. 
    8891073                  END IF 
     
    8911075            END DO 
    8921076            ! 
    893          END DO    
     1077         END DO 
    8941078         !! 
    8951079         !! BL89 search for ldwn 
    896          !! ----------------------------------------------------------           
    897          DO jk=2,jpka-1          
     1080         !! ---------------------------------------------------------- 
     1081         DO jk=2,jpka-1 
    8981082            ! 
    8991083            DO ji = 1, jpi 
     
    9011085               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a ) 
    9021086               ln_foundl(ji ) = .false. 
    903             END DO   
    904             !    
     1087            END DO 
     1088            ! 
    9051089            DO jkdwn=jk-1,1,-1 
    906                DO ji = 1, jpi              
     1090               DO ji = 1, jpi 
     1091                  zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 
     1092                  zbuoy2 = MAX( zbn2(ji,jj,jkdwn  ), rsmall ) 
    9071093                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  & 
    908                      &               * ( zbn2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 
    909                                        + zbn2(ji,jj,jkdwn  )*(ghw_abl(jk)-ghw_abl(jkdwn  )) )    
    910                   IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN  
     1094                     &               * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 
     1095                                       + zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn  )) ) 
     1096                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN 
    9111097                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 
    912                      zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )               
     1098                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  ) 
    9131099                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   & 
    914                         &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) )  
    915                      zldw(ji,jk) = zcff  
    916                      ln_foundl(ji) = .true.                
    917                   END IF                                    
    918                END DO 
    919             END DO 
    920             !       
     1100                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) ) 
     1101                     zldw(ji,jk) = zcff 
     1102                     zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  ) 
     1103                     ln_foundl(ji) = .true. 
     1104                  END IF 
     1105               END DO 
     1106            END DO 
     1107            ! 
    9211108         END DO 
    9221109 
    9231110         DO jk = 1, jpka 
    924             DO ji = 1, jpi   
    925                mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min )   
    926             END DO 
    927          END DO   
     1111            DO ji = 1, jpi 
     1112               !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp) 
     1113               zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     1114               mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 
     1115               mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min ) 
     1116            END DO 
     1117         END DO 
    9281118 
    9291119      END DO 
    930 #   undef zlup       
    931 #   undef zldw            
    932          ! 
    933       END SELECT       
     1120#   undef zlup 
     1121#   undef zldw 
     1122         ! 
     1123     CASE ( 3 )           ! Bougeault & Lacarrere 89 length-scale 
     1124         ! 
     1125#   define zlup zRH 
     1126#   define zldw zFC 
     1127! zCF is used for matrix inversion 
     1128! 
     1129       DO jj = 1, jpj      ! outer loop 
     1130          ! 
     1131          DO jk = 2, jpkam1 
     1132             DO ji = 1,jpi 
     1133                            zcff        = 1.0_wp / e3w_abl( jk )**2 
     1134                zdU         = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2 
     1135                zdV         = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2 
     1136                zsh2(ji,jk) = SQRT(zdU+zdV)   !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 
     1137                         ENDDO 
     1138                  ENDDO 
     1139          zsh2(:,      1) = zsh2( :,      2) 
     1140          zsh2(:,   jpka) = zsh2( :, jpkam1) 
     1141 
     1142                 DO ji = 1, jpi 
     1143            zcff              = zrough(ji,jj) * rn_Lsfc 
     1144                        zlup( ji,    1 )  = zcff 
     1145            zldw( ji,    1 )  = zcff 
     1146            zlup( ji, jpka ) = mxl_min 
     1147            zldw( ji, jpka ) = mxl_min 
     1148         END DO 
     1149 
     1150         DO jk = 2,jpka-1 
     1151            DO ji = 1, jpi 
     1152               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk) 
     1153               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1) 
     1154            END DO 
     1155         END DO 
     1156         !! 
     1157         !! BL89 search for lup 
     1158         !! ---------------------------------------------------------- 
     1159         DO jk=2,jpka-1 
     1160            ! 
     1161            DO ji = 1, jpi 
     1162               zCF(ji,1:jpka) = 0._wp 
     1163               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a ) 
     1164               ln_foundl(ji ) = .false. 
     1165            END DO 
     1166            ! 
     1167            DO jkup=jk+1,jpka-1 
     1168               DO ji = 1, jpi 
     1169                  zbuoy1             = MAX( zbn2(ji,jj,jkup  ), rsmall ) 
     1170                  zbuoy2             = MAX( zbn2(ji,jj,jkup-1), rsmall ) 
     1171                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) *                 & 
     1172                     &               ( zbuoy1*(ghw_abl(jkup  )-ghw_abl(jk))      & 
     1173                     &               + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) )    & 
     1174                                         &                            + 0.5_wp * e3t_abl(jkup) * rn_Rod *        & 
     1175                                         &               ( SQRT(tke_abl( ji, jj, jkup  , nt_a ))*zsh2(ji,jkup  ) & 
     1176                                         &               + SQRT(tke_abl( ji, jj, jkup-1, nt_a ))*zsh2(ji,jkup-1) ) 
     1177 
     1178                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 
     1179                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk) 
     1180                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 
     1181                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   & 
     1182                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) ) 
     1183                     zlup(ji,jk) = zcff 
     1184                     zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk) 
     1185                     ln_foundl(ji) = .true. 
     1186                  END IF 
     1187               END DO 
     1188            END DO 
     1189            ! 
     1190         END DO 
     1191         !! 
     1192         !! BL89 search for ldwn 
     1193         !! ---------------------------------------------------------- 
     1194         DO jk=2,jpka-1 
     1195            ! 
     1196            DO ji = 1, jpi 
     1197               zCF(ji,1:jpka) = 0._wp 
     1198               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a ) 
     1199               ln_foundl(ji ) = .false. 
     1200            END DO 
     1201            ! 
     1202            DO jkdwn=jk-1,1,-1 
     1203               DO ji = 1, jpi 
     1204                  zbuoy1             = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 
     1205                  zbuoy2             = MAX( zbn2(ji,jj,jkdwn  ), rsmall ) 
     1206                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  & 
     1207                     &               * (zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1))    & 
     1208                     &                 +zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn  )) )  & 
     1209                                         &                            + 0.5_wp * e3t_abl(jkup) * rn_Rod *          & 
     1210                                         &               ( SQRT(tke_abl( ji, jj, jkdwn+1, nt_a ))*zsh2(ji,jkdwn+1) & 
     1211                                         &               + SQRT(tke_abl( ji, jj, jkdwn  , nt_a ))*zsh2(ji,jkdwn  ) ) 
     1212 
     1213                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN 
     1214                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 
     1215                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  ) 
     1216                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   & 
     1217                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) ) 
     1218                     zldw(ji,jk) = zcff 
     1219                     zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  ) 
     1220                     ln_foundl(ji) = .true. 
     1221                  END IF 
     1222               END DO 
     1223            END DO 
     1224            ! 
     1225         END DO 
     1226 
     1227         DO jk = 1, jpka 
     1228            DO ji = 1, jpi 
     1229               !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp) 
     1230               zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     1231               mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 
     1232               mxld_abl( ji, jj, jk ) = MAX(  MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min ) 
     1233            END DO 
     1234         END DO 
     1235 
     1236      END DO 
     1237#   undef zlup 
     1238#   undef zldw 
     1239         ! 
     1240         ! 
     1241      END SELECT 
    9341242      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    9351243      !                            !  Finalize the computation of turbulent visc./diff. 
    9361244      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    937        
     1245 
    9381246      !------------- 
    9391247      DO jj = 1, jpj     ! outer loop 
    9401248      !------------- 
    941          DO jk = 1, jpka    
     1249         DO jk = 1, jpka 
    9421250            DO ji = 1, jpi  ! vector opt. 
    943                zcff              = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk )  & 
    944                   &                                       * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) )  
    945                zcff2             =  1. / ( 1. + zcff )   !<-- phi_z(z) 
    946                zcff              = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) ) 
    947 !!FL: MAX function probably useless because of the definition of mxl_min              
     1251               zcff  = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk )  & 
     1252               &     * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) ) 
     1253               zcff2 =  1. / ( 1. + zcff )   !<-- phi_z(z) 
     1254               zcff  = mxlm_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) ) 
     1255               !!FL: MAX function probably useless because of the definition of mxl_min 
    9481256               Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff         , avm_bak   ) 
    949                Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   )                               
    950             END DO 
    951          END DO 
    952       !-------------          
    953       END DO       
     1257               Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   ) 
     1258            END DO 
     1259         END DO 
     1260      !------------- 
     1261      END DO 
    9541262      !------------- 
    9551263 
     
    9691277      !! 
    9701278      !! --------------------------------------------------------------------- 
    971      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: msk    
    972      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d 
     1279      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: msk 
     1280      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d 
    9731281      INTEGER                                     :: ji,jj 
    974      REAL(wp)                                    :: smth_a, smth_b 
    975      REAL(wp), DIMENSION(jpi,jpj)                :: zdX,zdY,zFX,zFY 
    976      REAL(wp)                                    :: zumsk,zvmsk 
     1282      REAL(wp)                                    :: smth_a, smth_b 
     1283      REAL(wp), DIMENSION(jpi,jpj)                :: zdX,zdY,zFX,zFY 
     1284      REAL(wp)                                    :: zumsk,zvmsk 
    9771285      !! 
    9781286      !!========================================================= 
     
    9861294         zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji  ,jj ) ) * zumsk 
    9871295      END_2D 
    988        
    989      DO_2D_10_11 
     1296 
     1297      DO_2D_10_11 
    9901298         zvmsk = msk(ji,jj) * msk(ji,jj+1) 
    9911299         zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji  ,jj ) ) * zvmsk 
    992      END_2D 
    993        
    994      DO_2D_10_00 
     1300      END_2D 
     1301 
     1302      DO_2D_10_00 
    9951303         zFY ( ji, jj  ) =   zdY ( ji, jj   )                        & 
    9961304            & +  smth_a*  ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 ))   & 
    9971305            &            -  (zdX ( ji, jj   ) - zdX( ji-1, jj   ))  ) 
    998      END_2D 
     1306      END_2D 
    9991307 
    10001308      DO_2D_00_10 
     
    10101318  &                 +zFY( ji, jj ) - zFY( ji, jj-1 )  ) 
    10111319      END_2D 
    1012      !! 
     1320 
    10131321!--------------------------------------------------------------------------------------------------- 
    10141322   END SUBROUTINE smooth_pblh 
  • NEMO/trunk/src/ABL/ablrst.F90

    r12649 r13214  
    109109      CALL iom_delay_rst( 'WRITE', 'ABL', numraw )   ! save only abl delayed global communication variables 
    110110 
    111       ! Prognostic variables 
     111      ! Prognostic (after timestep + swap time indices = now timestep) variables 
    112112      CALL iom_rstput( iter, nitrst, numraw,   'u_abl',   u_abl(:,:,:,nt_n      ) ) 
    113113      CALL iom_rstput( iter, nitrst, numraw,   'v_abl',   v_abl(:,:,:,nt_n      ) ) 
     
    117117      CALL iom_rstput( iter, nitrst, numraw, 'avm_abl', avm_abl(:,:,:           ) ) 
    118118      CALL iom_rstput( iter, nitrst, numraw, 'avt_abl', avt_abl(:,:,:           ) ) 
    119       CALL iom_rstput( iter, nitrst, numraw, 'mxl_abl', mxl_abl(:,:,:           ) ) 
     119      CALL iom_rstput( iter, nitrst, numraw,'mxld_abl',mxld_abl(:,:,:           ) ) 
    120120      CALL iom_rstput( iter, nitrst, numraw,    'pblh',    pblh(:,:             ) ) 
    121121      ! 
     
    172172      CALL iom_get( numrar, jpdom_autoglo, 'avm_abl', avm_abl(:,:,:           ) ) 
    173173      CALL iom_get( numrar, jpdom_autoglo, 'avt_abl', avt_abl(:,:,:           ) ) 
    174       CALL iom_get( numrar, jpdom_autoglo, 'mxl_abl', mxl_abl(:,:,:           ) ) 
     174      CALL iom_get( numrar, jpdom_autoglo,'mxld_abl',mxld_abl(:,:,:           ) ) 
    175175      CALL iom_get( numrar, jpdom_autoglo,    'pblh',    pblh(:,:             ) ) 
    176176      CALL iom_delay_rst( 'READ', 'ABL', numrar )   ! read only abl delayed global communication variables 
  • NEMO/trunk/src/ABL/par_abl.F90

    r12814 r13214  
    2828   LOGICAL , PUBLIC            ::   ln_hpgls_frc   !: forcing of ABL winds by large-scale pressure gradient  
    2929   LOGICAL , PUBLIC            ::   ln_smth_pblh   !: smoothing of atmospheric PBL height  
     30   !LOGICAL , PUBLIC            ::   ln_topbc_neumann = .FALSE.  !: idealised testcases only 
    3031 
    31    LOGICAL           , PUBLIC ::   ln_rstart_abl    !: (de)activate abl restart 
    32    CHARACTER(len=256), PUBLIC ::   cn_ablrst_in     !: suffix of abl restart name (input) 
    33    CHARACTER(len=256), PUBLIC ::   cn_ablrst_out    !: suffix of abl restart name (output) 
    34    CHARACTER(len=256), PUBLIC ::   cn_ablrst_indir  !: abl restart input directory 
    35    CHARACTER(len=256), PUBLIC ::   cn_ablrst_outdir !: abl restart output directory 
     32   LOGICAL           , PUBLIC  ::   ln_rstart_abl    !: (de)activate abl restart 
     33   CHARACTER(len=256), PUBLIC  ::   cn_ablrst_in     !: suffix of abl restart name (input) 
     34   CHARACTER(len=256), PUBLIC  ::   cn_ablrst_out    !: suffix of abl restart name (output) 
     35   CHARACTER(len=256), PUBLIC  ::   cn_ablrst_indir  !: abl restart input directory 
     36   CHARACTER(len=256), PUBLIC  ::   cn_ablrst_outdir !: abl restart output directory 
    3637 
    3738   !!--------------------------------------------------------------------- 
     
    4647   REAL(wp), PUBLIC, PARAMETER ::   rn_Cek    = 258._wp                   !: Ekman constant for Richardson number  
    4748   REAL(wp), PUBLIC, PARAMETER ::   rn_epssfc = 1._wp / ( 1._wp + 2.8_wp * 2.8_wp ) 
    48    REAL(wp), PUBLIC            ::   rn_ceps                       !: namelist parameter 
    49    REAL(wp), PUBLIC            ::   rn_cm                         !: namelist parameter 
    50    REAL(wp), PUBLIC            ::   rn_ct                         !: namelist parameter 
    51    REAL(wp), PUBLIC            ::   rn_ce                         !: namelist parameter  
     49   REAL(wp), PUBLIC            ::   rn_Ceps                       !: namelist parameter 
     50   REAL(wp), PUBLIC            ::   rn_Cm                         !: namelist parameter 
     51   REAL(wp), PUBLIC            ::   rn_Ct                         !: namelist parameter 
     52   REAL(wp), PUBLIC            ::   rn_Ce                         !: namelist parameter  
    5253   REAL(wp), PUBLIC            ::   rn_Rod                        !: namelist parameter    
    5354   REAL(wp), PUBLIC            ::   rn_Sch     
     55   REAL(wp), PUBLIC            ::   rn_Esfc 
     56   REAL(wp), PUBLIC            ::   rn_Lsfc 
    5457   REAL(wp), PUBLIC            ::   mxl_min     
    5558   REAL(wp), PUBLIC            ::   rn_ldyn_min                   !: namelist parameter 
  • NEMO/trunk/src/ABL/sbcabl.F90

    r13208 r13214  
    7171         &                 ln_hpgls_frc, ln_geos_winds, nn_dyn_restore,           & 
    7272         &                 rn_ldyn_min , rn_ldyn_max, rn_ltra_min, rn_ltra_max,   & 
    73          &                 nn_amxl, rn_cm, rn_ct, rn_ce, rn_ceps, rn_Rod, rn_Ric, & 
     73         &                 nn_amxl, rn_Cm, rn_Ct, rn_Ce, rn_Ceps, rn_Rod, rn_Ric, & 
    7474         &                 ln_smth_pblh 
    7575      !!--------------------------------------------------------------------- 
    7676 
    77       ! Namelist namsbc_abl in reference namelist : ABL parameters 
     77                                        ! Namelist namsbc_abl in reference namelist : ABL parameters 
    7878      READ  ( numnam_ref, namsbc_abl, IOSTAT = ios, ERR = 901 ) 
    7979901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in reference namelist' ) 
    80       ! Namelist namsbc_abl in configuration namelist : ABL parameters 
     80      ! 
     81                                        ! Namelist namsbc_abl in configuration namelist : ABL parameters 
    8182      READ  ( numnam_cfg, namsbc_abl, IOSTAT = ios, ERR = 902 ) 
    8283902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in configuration namelist' ) 
     
    165166      rn_Sch  = rn_ce / rn_cm 
    166167      mxl_min = (avm_bak / rn_cm) / sqrt( tke_min ) 
     168      rn_Esfc =  1._wp / SQRT(rn_cm*rn_ceps) 
     169      rn_Lsfc = vkarmn * SQRT(SQRT(rn_cm*rn_ceps)) / rn_cm 
    167170 
    168171      IF(lwp) THEN 
     
    171174         WRITE(numout,*) '    ~~~~~~~~~~~' 
    172175         IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale ' 
    173          IF(nn_amxl==1) WRITE(numout,*) 'length-scale based on the distance to the PBL height ' 
     176         IF(nn_amxl==1) WRITE(numout,*) 'Modified Deardorff 80 length-scale ' 
     177         IF(nn_amxl==2) WRITE(numout,*) 'Bougeault and Lacarrere length-scale '       
     178         IF(nn_amxl==3) WRITE(numout,*) 'Rodier et al. length-scale '    
    174179         WRITE(numout,*) ' Minimum value of atmospheric TKE           = ',tke_min,' m^2 s^-2' 
    175180         WRITE(numout,*) ' Minimum value of atmospheric mixing length = ',mxl_min,' m' 
     
    178183         WRITE(numout,*) ' Constant for Schmidt number                = ',rn_Sch 
    179184         WRITE(numout,*) ' Constant for TKE dissipation               = ',rn_Ceps 
     185         WRITE(numout,*) ' Constant for TKE sfc boundary condition    = ',rn_Esfc 
     186         WRITE(numout,*) ' Constant for mxl sfc boundary condition    = ',rn_Lsfc 
    180187      END IF 
    181188 
     
    202209      ! ABL timestep 
    203210      rDt_abl = nn_fsbc * rn_Dt 
     211      IF(lwp) WRITE(numout,*) ' ABL timestep = ', rDt_abl,' s' 
    204212 
    205213      ! Check parameters for dynamics 
     
    248256         zcff         = 2._wp * omega * SIN( rad * 90._wp )   !++ fmax 
    249257         rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff) / zcff ) )**8 
    250          !!GS: alternative shape 
    251          !rest_eq(:,:) = SIN( 0.5_wp*rpi*(zcff - ABS(ff_t(:,:))) / (zcff - 3.e-5) )**8 
    252          !WHERE(ABS(ff_t(:,:)).LE.3.e-5) rest_eq(:,:) = 1._wp 
    253258      ELSE 
    254259         rest_eq(:,:) = 1._wp 
     
    271276         CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step 
    272277 
    273          u_abl(:,:,:,nt_n      ) = sf(jp_wndi)%fnow(:,:,:) 
    274          v_abl(:,:,:,nt_n      ) = sf(jp_wndj)%fnow(:,:,:) 
     278          u_abl(:,:,:,nt_n      ) = sf(jp_wndi)%fnow(:,:,:) 
     279          v_abl(:,:,:,nt_n      ) = sf(jp_wndj)%fnow(:,:,:) 
    275280         tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:) 
    276281         tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:) 
     
    279284         avm_abl(:,:,:          ) = avm_bak 
    280285         avt_abl(:,:,:          ) = avt_bak 
    281          mxl_abl(:,:,:          ) = mxl_min 
    282286         pblh   (:,:            ) = ghw_abl( 3 )  !<-- assume that the pbl contains 3 grid points 
    283287         u_abl  (:,:,:,nt_a     ) = 0._wp 
     
    285289         tq_abl (:,:,:,nt_a,:   ) = 0._wp 
    286290         tke_abl(:,:,:,nt_a     ) = 0._wp 
     291 
     292         mxlm_abl(:,:,:         ) = mxl_min 
     293         mxld_abl(:,:,:         ) = mxl_min 
    287294      ENDIF 
    288295 
  • NEMO/trunk/src/OCE/IOM/iom.F90

    r12649 r13214  
    11931193                           &         'we accept this case, even if there is a possible mix-up between z and time dimension' )    
    11941194                     idmspc = idmspc - 1 
    1195                   ELSE 
    1196                      CALL ctl_stop( TRIM(clinfo), 'To keep iom lisibility, when reading a '//clrankpv//'D array,'         ,   & 
    1197                         &                         'we do not accept data with '//cldmspc//' spatial dimensions',   & 
    1198                         &                         'Use ncwa -a to suppress the unnecessary dimensions' ) 
     1195                  !!GS: possibility to read 3D ABL atmopsheric forcing and use 1st level to force BULK simulation 
     1196                  !ELSE 
     1197                  !   CALL ctl_stop( TRIM(clinfo), 'To keep iom lisibility, when reading a '//clrankpv//'D array,',   & 
     1198                  !      &                         'we do not accept data with '//cldmspc//' spatial dimensions'  ,   & 
     1199                  !      &                         'Use ncwa -a to suppress the unnecessary dimensions' ) 
    11991200                  ENDIF 
    12001201            ENDIF 
  • NEMO/trunk/src/OCE/SBC/sbcblk.F90

    r13208 r13214  
    112112   REAL(wp)         ::   rn_zu     ! z(u)   : height of wind measurements 
    113113   ! 
    114    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice 
    115    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme) 
    116    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 
     114   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme and ABL) 
     115   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice 
     116   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 
    117117 
    118118   LOGICAL  ::   ln_skin_cs     ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB 
     
    121121   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 
    122122   LOGICAL  ::   ln_humi_rlh    ! humidity read in files ("sn_humi") is relative humidity     [%] if .true. !LB 
     123   LOGICAL  ::   ln_tpot        !!GS: flag to compute or not potential temperature 
    123124   ! 
    124125   INTEGER  ::   nhumi          ! choice of the bulk algorithm 
     
    181182         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm 
    182183         &                 cn_dir , rn_zqt, rn_zu,                                    & 
    183          &                 rn_pfac, rn_efac, ln_Cd_L12, ln_Cd_L15,                    & 
     184         &                 rn_pfac, rn_efac, ln_Cd_L12, ln_Cd_L15, ln_tpot,           & 
    184185         &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                          &   ! current feedback 
    185186         &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh  ! cool-skin / warm-layer !LB 
     
    626627         !#LB: because AGRIF hates functions that return something else than a scalar, need to 
    627628         !     use scalar version of gamma_moist() ... 
    628          DO_2D_11_11 
    629             ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 
    630          END_2D 
     629         IF( ln_tpot ) THEN 
     630            DO_2D_11_11 
     631               ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 
     632            END_2D 
     633         ELSE 
     634            ztpot = ptair(:,:) 
     635         ENDIF 
    631636      ENDIF 
    632637 
     
    970975         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   & 
    971976            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' ) 
    972       ELSE 
     977      ELSE ! ln_abl 
    973978         zztmp1 = 11637800.0_wp 
    974979         zztmp2 =    -5897.8_wp 
  • NEMO/trunk/tests/CANAL/MY_SRC/stpctl.F90

    r13136 r13214  
    1919   USE dom_oce         ! ocean space and time domain variables  
    2020   USE c1d             ! 1D vertical configuration 
    21    USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables 
    22    USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy 
    23    !   
    2421   USE diawri          ! Standard run outputs       (dia_wri_state routine) 
     22   ! 
    2523   USE in_out_manager  ! I/O manager 
    2624   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2725   USE lib_mpp         ! distributed memory computing 
    28    ! 
     26   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables 
     27   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy 
     28 
    2929   USE netcdf          ! NetCDF library 
    3030   IMPLICIT NONE 
     
    3333   PUBLIC stp_ctl           ! routine called by step.F90 
    3434 
    35    INTEGER                ::   nrunid   ! netcdf file id 
    36    INTEGER, DIMENSION(8)  ::   nvarid   ! netcdf variable id 
     35   INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus 
     36   LOGICAL  ::   lsomeoce 
    3737   !!---------------------------------------------------------------------- 
    3838   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    4242CONTAINS 
    4343 
    44    SUBROUTINE stp_ctl( kt, Kmm ) 
     44   SUBROUTINE stp_ctl( kt, Kbb, Kmm, kindic ) 
    4545      !!---------------------------------------------------------------------- 
    4646      !!                    ***  ROUTINE stp_ctl  *** 
     
    5050      !! ** Method  : - Save the time step in numstp 
    5151      !!              - Print it each 50 time steps 
    52       !!              - Stop the run IF problem encountered by setting nstop > 0 
     52      !!              - Stop the run IF problem encountered by setting indic=-3 
    5353      !!                Problems checked: |ssh| maximum larger than 10 m 
    5454      !!                                  |U|   maximum larger than 10 m/s  
     
    5757      !! ** Actions :   "time.step" file = last ocean time-step 
    5858      !!                "run.stat"  file = run statistics 
    59       !!                 nstop indicator sheared among all local domain 
     59      !!                nstop indicator sheared among all local domain (lk_mpp=T) 
    6060      !!---------------------------------------------------------------------- 
    6161      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index 
    62       INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index 
     62      INTEGER, INTENT(in   ) ::   Kbb, Kmm      ! ocean time level index 
     63      INTEGER, INTENT(inout) ::   kindic   ! error indicator 
    6364      !! 
    64       INTEGER                         ::   ji                                    ! dummy loop indices 
    65       INTEGER                         ::   idtime, istatus 
    66       INTEGER , DIMENSION(9)          ::   iareasum, iareamin, iareamax 
    67       INTEGER , DIMENSION(3,4)        ::   iloc                                  ! min/max loc indices 
    68       REAL(wp)                        ::   zzz                                   ! local real  
    69       REAL(wp), DIMENSION(9)          ::   zmax, zmaxlocal 
    70       LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns 
    71       LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk 
    72       CHARACTER(len=20)               ::   clname 
     65      INTEGER                ::   ji, jj, jk          ! dummy loop indices 
     66      INTEGER, DIMENSION(2)  ::   ih                  ! min/max loc indices 
     67      INTEGER, DIMENSION(3)  ::   iu, is1, is2        ! min/max loc indices 
     68      REAL(wp)               ::   zzz                 ! local real  
     69      REAL(wp), DIMENSION(9) ::   zmax 
     70      LOGICAL                ::   ll_wrtstp, ll_colruns, ll_wrtruns 
     71      CHARACTER(len=20) :: clname 
    7372      !!---------------------------------------------------------------------- 
    74       IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid 
    75       ! 
    76       ll_wrtstp  = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 
    77       ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1  
    78       ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm 
    79       ! 
    80       IF( kt == nit000 ) THEN 
    81          ! 
    82          IF( lwp ) THEN 
    83             WRITE(numout,*) 
    84             WRITE(numout,*) 'stp_ctl : time-stepping control' 
    85             WRITE(numout,*) '~~~~~~~' 
    86          ENDIF 
    87          !                                ! open time.step    ascii file, done only by 1st subdomain 
    88          IF( lwm )   CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    89          ! 
    90          IF( ll_wrtruns ) THEN 
    91             !                             ! open run.stat     ascii file, done only by 1st subdomain 
     73      ! 
     74      ll_wrtstp  = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 
     75      ll_colruns = ll_wrtstp .AND. ( sn_cfctl%l_runstat ) 
     76      ll_wrtruns = ll_colruns .AND. lwm 
     77      IF( kt == nit000 .AND. lwp ) THEN 
     78         WRITE(numout,*) 
     79         WRITE(numout,*) 'stp_ctl : time-stepping control' 
     80         WRITE(numout,*) '~~~~~~~' 
     81         !                                ! open time.step file 
     82         IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
     83         !                                ! open run.stat file(s) at start whatever 
     84         !                                ! the value of sn_cfctl%ptimincr 
     85         IF( lwm .AND. ( sn_cfctl%l_runstat ) ) THEN 
    9286            CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    93             !                             ! open run.stat.nc netcdf file, done only by 1st subdomain 
    9487            clname = 'run.stat.nc' 
    9588            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 
    96             istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, nrunid ) 
    97             istatus = NF90_DEF_DIM( nrunid, 'time', NF90_UNLIMITED, idtime ) 
    98             istatus = NF90_DEF_VAR( nrunid, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), nvarid(1) ) 
    99             istatus = NF90_DEF_VAR( nrunid,   'abs_u_max', NF90_DOUBLE, (/ idtime /), nvarid(2) ) 
    100             istatus = NF90_DEF_VAR( nrunid,       's_min', NF90_DOUBLE, (/ idtime /), nvarid(3) ) 
    101             istatus = NF90_DEF_VAR( nrunid,       's_max', NF90_DOUBLE, (/ idtime /), nvarid(4) ) 
    102             istatus = NF90_DEF_VAR( nrunid,       't_min', NF90_DOUBLE, (/ idtime /), nvarid(5) ) 
    103             istatus = NF90_DEF_VAR( nrunid,       't_max', NF90_DOUBLE, (/ idtime /), nvarid(6) ) 
     89            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, idrun ) 
     90            istatus = NF90_DEF_DIM( idrun, 'time', NF90_UNLIMITED, idtime ) 
     91            istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh ) 
     92            istatus = NF90_DEF_VAR( idrun,   'abs_u_max', NF90_DOUBLE, (/ idtime /), idu  ) 
     93            istatus = NF90_DEF_VAR( idrun,       's_min', NF90_DOUBLE, (/ idtime /), ids1 ) 
     94            istatus = NF90_DEF_VAR( idrun,       's_max', NF90_DOUBLE, (/ idtime /), ids2 ) 
     95            istatus = NF90_DEF_VAR( idrun,       't_min', NF90_DOUBLE, (/ idtime /), idt1 ) 
     96            istatus = NF90_DEF_VAR( idrun,       't_max', NF90_DOUBLE, (/ idtime /), idt2 ) 
    10497            IF( ln_zad_Aimp ) THEN 
    105                istatus = NF90_DEF_VAR( nrunid,   'Cf_max', NF90_DOUBLE, (/ idtime /), nvarid(7) ) 
    106                istatus = NF90_DEF_VAR( nrunid,'abs_wi_max',NF90_DOUBLE, (/ idtime /), nvarid(8) ) 
     98               istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1 ) 
     99               istatus = NF90_DEF_VAR( idrun,       'Cf_max', NF90_DOUBLE, (/ idtime /), idc1 ) 
    107100            ENDIF 
    108             istatus = NF90_ENDDEF(nrunid) 
    109          ENDIF 
    110          !     
    111       ENDIF 
    112       ! 
    113       !                                   !==              write current time step              ==! 
    114       !                                   !==  done only by 1st subdomain at writting timestep  ==! 
    115       IF( lwm .AND. ll_wrtstp ) THEN 
     101            istatus = NF90_ENDDEF(idrun) 
     102            zmax(8:9) = 0._wp    ! initialise to zero in case ln_zad_Aimp option is not in use 
     103         ENDIF 
     104      ENDIF 
     105      IF( kt == nit000 )   lsomeoce = COUNT( ssmask(:,:) == 1._wp ) > 0 
     106      ! 
     107      IF(lwm .AND. ll_wrtstp) THEN        !==  current time step  ==!   ("time.step" file) 
    116108         WRITE ( numstp, '(1x, i8)' )   kt 
    117109         REWIND( numstp ) 
    118110      ENDIF 
    119       !                                   !==            test of local extrema           ==! 
    120       !                                   !==  done by all processes at every time step  ==! 
    121       ! 
    122       ! define zmax default value. needed for land processors 
    123       IF( ll_colruns ) THEN    ! default value: must not be kept when calling mpp_max -> must be as small as possible 
    124          zmax(:) = -HUGE(1._wp) 
    125       ELSE                     ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...) 
    126          zmax(:) =  0._wp 
    127          zmax(3) = -1._wp      ! avoid salinity minimum at 0. 
    128       ENDIF 
    129       ! 
    130       llmsk(:,:,1) = ssmask(:,:) == 1._wp 
    131       IF( COUNT( llmsk(:,:,1) ) > 0 ) THEN   ! avoid huge values sent back for land processors... 
    132          IF( ll_wd ) THEN 
    133             zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max 
    134          ELSE 
    135             zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm)           ), mask = llmsk(:,:,1) )   ! ssh max 
    136          ENDIF 
    137       ENDIF 
    138       zmax(2) = MAXVAL( ABS( uu(:,:,:,Kmm) ) )                                       ! velocity max (zonal only) 
    139       llmsk(:,:,:) = tmask(:,:,:) == 1._wp 
    140       IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors... 
    141          zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     ! minus salinity max 
    142          zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     !       salinity max 
    143          IF( ll_colruns .OR. jpnij == 1 ) THEN     ! following variables are used only in the netcdf file 
    144             zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  ! minus temperature max 
    145             zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  !       temperature max 
    146             IF( ln_zad_Aimp ) THEN 
    147                zmax(7) = MAXVAL(   Cu_adv(:,:,:)   , mask = llmsk )                  ! partitioning coeff. max 
    148                llmsk(:,:,:) = wmask(:,:,:) == 1._wp 
    149                IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors... 
    150                   zmax(8) = MAXVAL(ABS( wi(:,:,:) ), mask = llmsk )                  ! implicit vertical vel. max 
    151                ENDIF 
    152             ENDIF 
    153          ENDIF 
    154       ENDIF 
    155       zmax(9) = REAL( nstop, wp )                                              ! stop indicator 
    156       !                                   !==               get global extrema             ==! 
    157       !                                   !==  done by all processes if writting run.stat  ==! 
     111      ! 
     112      !                                   !==  test of extrema  ==! 
     113      IF( ll_wd ) THEN 
     114         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) )  )        ! ssh max  
     115      ELSE 
     116         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) )  )                               ! ssh max 
     117      ENDIF 
     118      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) )  )                                  ! velocity max (zonal only) 
     119      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max 
     120      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       salinity max 
     121      zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max 
     122      zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       temperature max 
     123      zmax(7) = REAL( nstop , wp )                                            ! stop indicator 
     124      IF( ln_zad_Aimp ) THEN 
     125         zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 
     126         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max 
     127      ENDIF 
     128      ! 
    158129      IF( ll_colruns ) THEN 
    159          zmaxlocal(:) = zmax(:) 
    160130         CALL mpp_max( "stpctl", zmax )          ! max over the global domain 
    161          nstop = NINT( zmax(9) )                 ! update nstop indicator (now sheared among all local domains) 
    162       ENDIF 
    163       !                                   !==              write "run.stat" files              ==! 
    164       !                                   !==  done only by 1st subdomain at writting timestep  ==! 
     131         nstop = NINT( zmax(7) )                 ! nstop indicator sheared among all local domains 
     132      ENDIF 
     133      !                                   !==  run statistics  ==!   ("run.stat" files) 
    165134      IF( ll_wrtruns ) THEN 
    166135         WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 
    167          istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) ) 
    168          istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) ) 
    169          istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) ) 
    170          istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) ) 
    171          istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) ) 
    172          istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) ) 
     136         istatus = NF90_PUT_VAR( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) ) 
     137         istatus = NF90_PUT_VAR( idrun,   idu, (/ zmax(2)/), (/kt/), (/1/) ) 
     138         istatus = NF90_PUT_VAR( idrun,  ids1, (/-zmax(3)/), (/kt/), (/1/) ) 
     139         istatus = NF90_PUT_VAR( idrun,  ids2, (/ zmax(4)/), (/kt/), (/1/) ) 
     140         istatus = NF90_PUT_VAR( idrun,  idt1, (/-zmax(5)/), (/kt/), (/1/) ) 
     141         istatus = NF90_PUT_VAR( idrun,  idt2, (/ zmax(6)/), (/kt/), (/1/) ) 
    173142         IF( ln_zad_Aimp ) THEN 
    174             istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) ) 
    175             istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) ) 
    176          ENDIF 
    177          IF( kt == nitend )   istatus = NF90_CLOSE(nrunid) 
     143            istatus = NF90_PUT_VAR( idrun,  idw1, (/ zmax(8)/), (/kt/), (/1/) ) 
     144            istatus = NF90_PUT_VAR( idrun,  idc1, (/ zmax(9)/), (/kt/), (/1/) ) 
     145         ENDIF 
     146         IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun) 
     147         IF( kt == nitend         ) istatus = NF90_CLOSE(idrun) 
    178148      END IF 
    179       !                                   !==               error handling               ==! 
    180       !                                   !==  done by all processes at every time step  ==! 
    181       ! 
    182       IF(   zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m ) 
    183          &  zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s) 
    184 !!$         &  zmax(3) >=   0._wp .OR.   &                   ! negative or zero sea surface salinity 
    185 !!$         &  zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 ) 
    186 !!$         &  zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice) 
    187          &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests 
    188          &  ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests 
     149      !                                   !==  error handling  ==! 
     150      IF( ( sn_cfctl%l_glochk .OR. lsomeoce ) .AND. (   &  ! domain contains some ocean points, check for sensible ranges 
     151         &  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m ) 
     152         &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
     153!!$         &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity 
     154!!$         &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 ) 
     155!!$         &  zmax(4) <    0._wp .OR.   &                    ! too large sea surface salinity (keep this line for sea-ice) 
     156         &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN   ! NaN encounter in the tests 
     157         IF( lk_mpp .AND. sn_cfctl%l_glochk ) THEN 
     158            ! have use mpp_max (because sn_cfctl%l_glochk=.T. and distributed) 
     159            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm))        , ssmask(:,:)  , zzz, ih  ) 
     160            CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm))          , umask (:,:,:), zzz, iu  ) 
     161            CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is1 ) 
     162            CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is2 ) 
     163         ELSE 
     164            ! find local min and max locations 
     165            ih(:)  = MAXLOC( ABS( ssh(:,:,Kmm)   )                              ) + (/ nimpp - 1, njmpp - 1    /) 
     166            iu(:)  = MAXLOC( ABS( uu  (:,:,:,Kmm) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     167            is1(:) = MINLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     168            is2(:) = MAXLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     169         ENDIF 
     170          
     171         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
     172         WRITE(ctmp2,9100) kt,   zmax(1), ih(1) , ih(2) 
     173         WRITE(ctmp3,9200) kt,   zmax(2), iu(1) , iu(2) , iu(3) 
     174         WRITE(ctmp4,9300) kt, - zmax(3), is1(1), is1(2), is1(3) 
     175         WRITE(ctmp5,9400) kt,   zmax(4), is2(1), is2(2), is2(3) 
     176         WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort.nc file' 
     177          
     178         CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file 
     179          
     180         IF( .NOT. sn_cfctl%l_glochk ) THEN 
     181            WRITE(ctmp8,*) 'E R R O R message from sub-domain: ', narea 
     182            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ctmp6 ) 
     183         ELSE 
     184            CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6, ' ' ) 
     185         ENDIF 
     186 
     187         kindic = -3 
    189188         ! 
    190          iloc(:,:) = 0 
    191          IF( ll_colruns ) THEN   ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc 
    192             ! first: close the netcdf file, so we can read it 
    193             IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(nrunid) 
    194             ! get global loc on the min/max 
    195             CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,         Kmm)), ssmask(:,:  ), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F  
    196             CALL mpp_maxloc( 'stpctl', ABS( uu(:,:,:,       Kmm)),  umask(:,:,:), zzz, iloc(1:3,2) ) 
    197             CALL mpp_minloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) ,  tmask(:,:,:), zzz, iloc(1:3,3) ) 
    198             CALL mpp_maxloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) ,  tmask(:,:,:), zzz, iloc(1:3,4) ) 
    199             ! find which subdomain has the max. 
    200             iareamin(:) = jpnij+1   ;   iareamax(:) = 0   ;   iareasum(:) = 0 
    201             DO ji = 1, 9 
    202                IF( zmaxlocal(ji) == zmax(ji) ) THEN 
    203                   iareamin(ji) = narea   ;   iareamax(ji) = narea   ;   iareasum(ji) = 1 
    204                ENDIF 
    205             END DO 
    206             CALL mpp_min( "stpctl", iareamin )         ! min over the global domain 
    207             CALL mpp_max( "stpctl", iareamax )         ! max over the global domain 
    208             CALL mpp_sum( "stpctl", iareasum )         ! sum over the global domain 
    209          ELSE                    ! find local min and max locations: 
    210             ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc 
    211             iloc(1:2,1) = MAXLOC( ABS( ssh(:,:,         Kmm)), mask = ssmask(:,:  ) == 1._wp ) + (/ nimpp - 1, njmpp - 1    /) 
    212             iloc(1:3,2) = MAXLOC( ABS(  uu(:,:,:,       Kmm)), mask =  umask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    213             iloc(1:3,3) = MINLOC(       ts(:,:,:,jp_sal,Kmm) , mask =  tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    214             iloc(1:3,4) = MAXLOC(       ts(:,:,:,jp_sal,Kmm) , mask =  tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    215             iareamin(:) = narea   ;   iareamax(:) = narea   ;   iareasum(:) = 1         ! this is local information 
    216          ENDIF 
    217          ! 
    218          WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
    219          CALL wrt_line( ctmp2, kt, '|ssh| max',  zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 
    220          CALL wrt_line( ctmp3, kt, '|U|   max',  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 
    221          CALL wrt_line( ctmp4, kt, 'Sal   min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 
    222          CALL wrt_line( ctmp5, kt, 'Sal   max',  zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 
    223          IF( Agrif_Root() ) THEN 
    224             WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files' 
    225          ELSE 
    226             WRITE(ctmp6,*) '      ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files' 
    227          ENDIF 
    228          ! 
    229          CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file 
    230          ! 
    231          IF( ll_colruns .or. jpnij == 1 ) THEN   ! all processes synchronized -> use lwp to print in opened ocean.output files 
    232             IF(lwp) THEN   ;   CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 ) 
    233             ELSE           ;   nstop = MAX(1, nstop)   ! make sure nstop > 0 (automatically done when calling ctl_stop) 
    234             ENDIF 
    235          ELSE                                    ! only mpi subdomains with errors are here -> STOP now 
    236             CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 ) 
    237          ENDIF 
    238          ! 
    239       ENDIF 
    240       ! 
    241       IF( nstop > 0 ) THEN                                                  ! an error was detected and we did not abort yet... 
    242          ngrdstop = Agrif_Fixed()                                           ! store which grid got this error 
    243          IF( .NOT. ll_colruns .AND. jpnij > 1 )   CALL ctl_stop( 'STOP' )   ! we must abort here to avoid MPI deadlock 
    244       ENDIF 
    245       ! 
     189      ENDIF 
     190      ! 
     1919100  FORMAT (' kt=',i8,'   |ssh| max: ',1pg11.4,', at  i j  : ',2i5) 
     1929200  FORMAT (' kt=',i8,'   |U|   max: ',1pg11.4,', at  i j k: ',3i5) 
     1939300  FORMAT (' kt=',i8,'   S     min: ',1pg11.4,', at  i j k: ',3i5) 
     1949400  FORMAT (' kt=',i8,'   S     max: ',1pg11.4,', at  i j k: ',3i5) 
    2461959500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 
    247196      ! 
    248197   END SUBROUTINE stp_ctl 
    249  
    250  
    251    SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax ) 
    252       !!---------------------------------------------------------------------- 
    253       !!                     ***  ROUTINE wrt_line  *** 
    254       !! 
    255       !! ** Purpose :   write information line 
    256       !! 
    257       !!---------------------------------------------------------------------- 
    258       CHARACTER(len=*),      INTENT(  out) ::   cdline 
    259       CHARACTER(len=*),      INTENT(in   ) ::   cdprefix 
    260       REAL(wp),              INTENT(in   ) ::   pval 
    261       INTEGER, DIMENSION(3), INTENT(in   ) ::   kloc 
    262       INTEGER,               INTENT(in   ) ::   kt, ksum, kmin, kmax 
    263       ! 
    264       CHARACTER(len=80) ::   clsuff 
    265       CHARACTER(len=9 ) ::   clkt, clsum, clmin, clmax 
    266       CHARACTER(len=9 ) ::   cli, clj, clk 
    267       CHARACTER(len=1 ) ::   clfmt 
    268       CHARACTER(len=4 ) ::   cl4   ! needed to be able to compile with Agrif, I don't know why 
    269       INTEGER           ::   ifmtk 
    270       !!---------------------------------------------------------------------- 
    271       WRITE(clkt , '(i9)') kt 
    272        
    273       WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij  ,wp))) + 1     ! how many digits to we need to write ? (we decide max = 9) 
    274       !!! WRITE(clsum, '(i'//clfmt//')') ksum                   ! this is creating a compilation error with AGRIF 
    275       cl4 = '(i'//clfmt//')'   ;   WRITE(clsum, cl4) ksum 
    276       WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1    ! how many digits to we need to write ? (we decide max = 9) 
    277       cl4 = '(i'//clfmt//')'   ;   WRITE(clmin, cl4) kmin-1 
    278                                    WRITE(clmax, cl4) kmax-1 
    279       ! 
    280       WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1      ! how many digits to we need to write jpiglo? (we decide max = 9) 
    281       cl4 = '(i'//clfmt//')'   ;   WRITE(cli, cl4) kloc(1)      ! this is ok with AGRIF 
    282       WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1      ! how many digits to we need to write jpjglo? (we decide max = 9) 
    283       cl4 = '(i'//clfmt//')'   ;   WRITE(clj, cl4) kloc(2)      ! this is ok with AGRIF 
    284       ! 
    285       IF( ksum == 1 ) THEN   ;   WRITE(clsuff,9100) TRIM(clmin) 
    286       ELSE                   ;   WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax) 
    287       ENDIF 
    288       IF(kloc(3) == 0) THEN 
    289          ifmtk = INT(LOG10(REAL(jpk,wp))) + 1                   ! how many digits to we need to write jpk? (we decide max = 9) 
    290          clk = REPEAT(' ', ifmtk)                               ! create the equivalent in blank string 
    291          WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff) 
    292       ELSE 
    293          WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1      ! how many digits to we need to write jpk? (we decide max = 9) 
    294          !!! WRITE(clk, '(i'//clfmt//')') kloc(3)               ! this is creating a compilation error with AGRIF 
    295          cl4 = '(i'//clfmt//')'   ;   WRITE(clk, cl4) kloc(3)   ! this is ok with AGRIF 
    296          WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj),    TRIM(clk), TRIM(clsuff) 
    297       ENDIF 
    298       ! 
    299 9100  FORMAT('MPI rank ', a) 
    300 9200  FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a) 
    301 9300  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j   ', a, ' ', a, ' ', a, ' ', a) 
    302 9400  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a) 
    303       ! 
    304    END SUBROUTINE wrt_line 
    305  
    306198 
    307199   !!====================================================================== 
  • NEMO/trunk/tests/STATION_ASF/EXPREF/launch_sasf.sh

    r13132 r13214  
    11#!/bin/bash 
    22 
    3 ################################################################ 
    4 # 
    5 # Script to launch a set of STATION_ASF simulations 
    6 # 
    7 # L. Brodeau, 2020 
    8 # 
    9 ################################################################ 
     3# NEMO directory where to fetch compiled STATION_ASF nemo.exe + setup: 
     4NEMO_DIR=`pwd | sed -e "s|/tests/STATION_ASF/EXPREF||g"` 
    105 
    11 # What directory inside "tests" actually contains the compiled "nemo.exe" for STATION_ASF ? 
     6echo "Using NEMO_DIR=${NEMO_DIR}" 
     7 
     8# what directory inside "tests" actually contains the compiled test-case? 
    129TC_DIR="STATION_ASF2" 
    1310 
    14 # DATA_IN_DIR => Directory containing sea-surface + atmospheric forcings 
     11# => so the executable to use is: 
     12NEMO_EXE="${NEMO_DIR}/tests/${TC_DIR}/BLD/bin/nemo.exe" 
     13 
     14# Directory where to run the simulation: 
     15WORK_DIR="${HOME}/tmp/STATION_ASF" 
     16 
     17 
     18# FORC_DIR => Directory containing sea-surface + atmospheric forcings 
    1519#             (get it there https://drive.google.com/file/d/1MxNvjhRHmMrL54y6RX7WIaM9-LGl--ZP/): 
    1620if [ `hostname` = "merlat"        ]; then 
    17     DATA_IN_DIR="/MEDIA/data/STATION_ASF/input_data_STATION_ASF_2016-2018" 
     21    FORC_DIR="/MEDIA/data/STATION_ASF/input_data_STATION_ASF_2016-2018" 
    1822elif [ `hostname` = "luitel"        ]; then 
    19     DATA_IN_DIR="/data/gcm_setup/STATION_ASF/input_data_STATION_ASF_2016-2018" 
     23    FORC_DIR="/data/gcm_setup/STATION_ASF/input_data_STATION_ASF_2016-2018" 
    2024elif [ `hostname` = "ige-meom-cal1" ]; then 
    21     DATA_IN_DIR="/mnt/meom/workdir/brodeau/STATION_ASF/input_data_STATION_ASF_2016-2018" 
     25    FORC_DIR="/mnt/meom/workdir/brodeau/STATION_ASF/input_data_STATION_ASF_2016-2018" 
    2226elif [ `hostname` = "salvelinus" ]; then 
    23     DATA_IN_DIR="/opt/data/STATION_ASF/input_data_STATION_ASF_2016-2018" 
     27    FORC_DIR="/opt/data/STATION_ASF/input_data_STATION_ASF_2016-2018" 
    2428else 
    25     echo "Oops! We don't know `hostname` yet! Define 'DATA_IN_DIR' in the script!"; exit  
     29    echo "Boo!"; exit 
    2630fi 
    27  
    28 expdir=`basename ${PWD}`; # we expect "EXPREF" or "EXP00" normally... 
    29  
    30 # NEMOGCM root directory where to fetch compiled STATION_ASF nemo.exe + setup: 
    31 NEMO_WRK_DIR=`pwd | sed -e "s|/tests/STATION_ASF/${expdir}||g"` 
    32  
    33 # Directory where to run the simulation: 
    34 PROD_DIR="${HOME}/tmp/STATION_ASF" 
     31#====================== 
     32mkdir -p ${WORK_DIR} 
    3533 
    3634 
    37 ####### End of normal user configurable section ####### 
     35if [ ! -f ${NEMO_EXE} ]; then echo " Mhhh, no compiled nemo.exe found into ${NEMO_DIR}/tests/STATION_ASF/BLD/bin !"; exit; fi 
    3836 
    39 #================================================================================ 
    40  
    41 # NEMO executable to use is: 
    42 NEMO_EXE="${NEMO_WRK_DIR}/tests/${TC_DIR}/BLD/bin/nemo.exe" 
    43  
    44  
    45 echo "###########################################################" 
    46 echo "#        S T A T I O N   A i r  -  S e a   F l u x        #" 
    47 echo "###########################################################" 
    48 echo 
    49 echo " We shall work in here: ${STATION_ASF_DIR}/" 
    50 echo " NEMOGCM   work    depository is: ${NEMO_WRK_DIR}/" 
    51 echo "   ==> NEMO EXE to use: ${NEMO_EXE}" 
    52 echo " Input forcing data into: ${DATA_IN_DIR}/" 
    53 echo " Production will be done into: ${PROD_DIR}/" 
    54 echo 
    55  
    56 mkdir -p ${PROD_DIR} 
    57  
    58 if [ ! -f ${NEMO_EXE} ]; then echo " Mhhh, no compiled 'nemo.exe' found into `dirname ${NEMO_EXE}` !"; exit; fi 
    59  
    60 echo 
    61 echo " *** Using the following NEMO executable:" 
    62 echo "  ${NEMO_EXE} " 
    63 echo 
    64  
    65 NEMO_EXPREF="${NEMO_WRK_DIR}/tests/STATION_ASF/EXPREF" 
     37NEMO_EXPREF="${NEMO_DIR}/tests/STATION_ASF/EXPREF" 
    6638if [ ! -d ${NEMO_EXPREF} ]; then echo " Mhhh, no EXPREF directory ${NEMO_EXPREF} !"; exit; fi 
    6739 
    68 rsync -avP ${NEMO_EXE}          ${PROD_DIR}/ 
     40rsync -avP ${NEMO_EXE}          ${WORK_DIR}/ 
    6941 
    7042for ff in "context_nemo.xml" "domain_def_nemo.xml" "field_def_nemo-oce.xml" "file_def_nemo-oce.xml" "grid_def_nemo.xml" "iodef.xml" "namelist_ref"; do 
    7143    if [ ! -f ${NEMO_EXPREF}/${ff} ]; then echo " Mhhh, ${ff} not found into ${NEMO_EXPREF} !"; exit; fi 
    72     rsync -avPL ${NEMO_EXPREF}/${ff} ${PROD_DIR}/ 
     44    rsync -avPL ${NEMO_EXPREF}/${ff} ${WORK_DIR}/ 
    7345done 
    7446 
    7547# Copy forcing to work directory: 
    76 rsync -avP ${DATA_IN_DIR}/Station_PAPA_50N-145W*.nc ${PROD_DIR}/ 
     48rsync -avP ${FORC_DIR}/Station_PAPA_50N-145W*.nc ${WORK_DIR}/ 
    7749 
    7850for CASE in "ECMWF" "COARE3p6" "NCAR" "ECMWF-noskin" "COARE3p6-noskin"; do 
     
    8658    scase=`echo "${CASE}" | tr '[:upper:]' '[:lower:]'` 
    8759 
    88     rm -f ${PROD_DIR}/namelist_cfg 
    89     rsync -avPL ${NEMO_EXPREF}/namelist_${scase}_cfg ${PROD_DIR}/namelist_cfg 
     60    rm -f ${WORK_DIR}/namelist_cfg 
     61    rsync -avPL ${NEMO_EXPREF}/namelist_${scase}_cfg ${WORK_DIR}/namelist_cfg 
    9062 
    91     cd ${PROD_DIR}/ 
     63    cd ${WORK_DIR}/ 
    9264    echo 
    9365    echo "Launching NEMO !" 
Note: See TracChangeset for help on using the changeset viewer.