Changeset 12588


Ignore:
Timestamp:
2020-03-23T18:21:59+01:00 (2 weeks ago)
Author:
gsamson
Message:

revised ABL model version including:

  • albmod cleaning
  • new abl mixing length option (nn_amxl = 3)
  • use rho_air function from aerobulk everywhere
  • remove mxl_abl (replaced by master (mxlm_abl) and dissipative (mxlm_abl) mixing lengths)
  • temporary flag "ln_tpot" to disable potential temperature computation in sbcblk
Location:
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/abl.F90

    r12489 r12588  
    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] 
     
    3839   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   msk_abl 
    3940   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   rest_eq 
     41    
     42   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::    cft_abl 
     43   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:  )     ::   taux_abl 
     44   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:  )     ::   tauy_abl   
    4045   ! 
    4146   INTEGER , PUBLIC :: nt_n, nt_a       !: now / after indices (equal 1 or 2) 
     
    5560      !!---------------------------------------------------------------------- 
    5661      ! 
    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 ) 
     62      ALLOCATE( u_abl   (1:jpi,1:jpj,1:jpka,jptime     ), & 
     63         &      v_abl   (1:jpi,1:jpj,1:jpka,jptime     ), & 
     64         &      tq_abl  (1:jpi,1:jpj,1:jpka,jptime,jptq), & 
     65         &      tke_abl (1:jpi,1:jpj,1:jpka,jptime     ), & 
     66         &      avm_abl (1:jpi,1:jpj,1:jpka            ), & 
     67         &      avt_abl (1:jpi,1:jpj,1:jpka            ), & 
     68         &      mxld_abl(1:jpi,1:jpj,1:jpka            ), &          
     69         &      mxlm_abl(1:jpi,1:jpj,1:jpka            ), &  
     70         &      cft_abl (1:jpi,1:jpj,1:jpka            ), &  
     71         &      fft_abl (1:jpi,1:jpj                   ), & 
     72         &      pblh    (1:jpi,1:jpj                   ), &          
     73         &      taux_abl(1:jpi,1:jpj                   ), &    
     74         &      tauy_abl(1:jpi,1:jpj                   ), &         
     75         &      msk_abl (1:jpi,1:jpj                   ), & 
     76         &      rest_eq (1:jpi,1:jpj                   ), &          
     77         &      e3t_abl (1:jpka), e3w_abl(1:jpka)       , & 
     78         &      ght_abl (1:jpka), ghw_abl(1:jpka)       , STAT=ierr ) 
    6979         ! 
    7080      abl_alloc = ierr 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablmod.F90

    r12489 r12588  
    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 
    2020   USE sbcblk         ! use rn_?fac 
     
    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" 
     
    3939!=================================================================================================== 
    4040   SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, &                            ! in 
    41               &            pu_dta, pv_dta, pt_dta, pq_dta,    &  
     41              &            pu_dta, pv_dta, pt_dta, pq_dta,    & 
    4242              &            pslp_dta, pgu_dta, pgv_dta,        & 
    4343              &            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            &   ) 
     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      
     101      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (V-point) 
     102#endif 
    103103      ! 
    104104      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    
     105      REAL(wp), DIMENSION(1:jpi,2:jpka  )        ::   zCF 
    107106      ! 
    108107      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_a 
     
    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 = .FALSE. 
     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(Cdn_oce(ji,jj)) ) !<-- 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 
     175               DO ji = 2, jpim1 
     176               !DO ji = 1,jpi  !!GS: to be checked 
    174177                  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 
     
    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             
     209            DO ji = 1,jpi 
    207210               zcff                       =  1._wp / z_elem_b( ji, 2 ) 
    208211               zCF   (ji,   2           ) = - zcff * z_elem_c( ji, 2 ) 
    209212               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             
     213            END DO 
     214 
     215            DO jk = 3, jpka 
     216               DO ji = 1,jpi 
    214217                  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 ) 
     
    218221               END DO 
    219222            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   
     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      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     240      IF( SemiImp_Cor ) THEN 
     241 
     242      !------------- 
    239243      DO jk = 2, jpka    ! outer loop 
    240244      !------------- 
    241          !              
     245         ! 
    242246         ! Advance u_abl & v_abl to time n+1 
    243247         DO_2D_11_11 
    244248            zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl )  ! (f dt)**2 
    245     
     249 
    246250            u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(  & 
    247251               &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n )    & 
    248252               &                 +  rDt_abl * fft_abl(ji, jj) * v_abl ( ji , jj  , jk, nt_n ) )  & 
    249253               &                               / (1._wp + gamma_Cor*gamma_Cor*zcff) 
    250                 
     254 
    251255            v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(  & 
    252256               &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n )   & 
    253257               &                 -  rDt_abl * fft_abl(ji, jj) * u_abl ( ji   , jj, jk, nt_n )  ) & 
    254                &                                / (1._wp + gamma_Cor*gamma_Cor*zcff)                 
     258               &                                / (1._wp + gamma_Cor*gamma_Cor*zcff) 
    255259         END_2D 
    256          !                                    
     260         ! 
    257261      !------------- 
    258262      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      !------------- 
     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      ELSE 
     279 
     280      !------------- 
     281      DO jk = 2, jpka    ! outer loop 
     282      !------------- 
     283         ! 
     284         IF( MOD( kt, 2 ) == 0 ) then 
     285                 ! Advance u_abl & v_abl to time n+1 
     286         DO jj = 1, jpj 
     287            DO ji = 1, jpi 
     288               zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_n ) - pgv_dta(ji  ,jj  ,jk)  ) 
     289               u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff 
     290               zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_a ) - pgu_dta(ji  ,jj  ,jk)  ) 
     291               v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff ) 
     292               u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  u_abl( ji, jj, jk, nt_a ) 
     293            END DO 
     294         END DO 
     295         ! 
     296         ELSE 
     297         ! 
     298         DO jj = 1, jpj 
     299            DO ji = 1, jpi 
     300               zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_n ) - pgu_dta(ji  ,jj  ,jk)  ) 
     301               v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff 
     302               zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_a ) - pgv_dta(ji  ,jj  ,jk)  ) 
     303               u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff ) 
     304               v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  v_abl( ji, jj, jk, nt_a ) 
     305            END DO 
     306         END DO 
     307         ! 
     308         ENDIF 
     309         ! 
     310      !------------- 
     311      END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 
     312      !------------- 
     313 
     314      ENDIF 
     315 
     316      !------------- 
     317      ! 
     318      IF( ln_hpgls_frc ) THEN 
     319         DO jj = 1, jpj    ! outer loop 
    263320            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  
     321               DO ji = 1, jpi 
    279322                  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)   
     323                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk) 
    281324               ENDDO 
    282325            ENDDO 
    283          ENDDO  
     326         ENDDO 
    284327      END IF 
    285       
    286       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    287       !                            !  4 *** Advance u,v to time n+1  
    288       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  
    289       ! 
    290       !  Vertical diffusion for u_abl      
     328 
     329      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     330      !                            !  4 *** Advance u,v to time n+1 
     331      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     332      ! 
     333      !  Vertical diffusion for u_abl 
    291334      !------------- 
    292335      DO jj = 1, jpj    ! outer loop 
    293       !-------------     
     336      !------------- 
    294337 
    295338         DO jk = 3, jpkam1 
    296             DO ji = 1, jpi   
     339            DO ji = 1, jpi 
    297340               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                 
     341               z_elem_c( ji,     jk ) = - rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal 
    299342               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                             !       diagonal 
    300343            END DO 
    301          END DO   
    302              
    303          DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi)             
     344         END DO 
     345 
     346         DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi) 
    304347            !++ Surface boundary condition 
    305348            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              
     349            z_elem_c( ji,     2    ) = - rDt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   ) 
     350            ! 
     351            zztmp1  = pcd_du(ji, jj) 
     352            zztmp2  = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) ) * rn_vfac 
     353#if defined key_si3 
    311354            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          
     355            zzice  = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) ) * rn_vfac 
     356            zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 
     357#endif 
     358            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rDt_abl * zztmp1 
    316359            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  
     360 
     361            IF( ln_topbc_neumann ) THEN 
     362               !++ Top Neumann B.C. 
     363               z_elem_a( ji,     jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
     364               z_elem_c( ji,     jpka ) = 0._wp 
     365               z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka ) 
     366               !u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * u_abl   ( ji, jj, jpka, nt_a ) 
     367            ELSE 
     368               !++ Top Dirichlet B.C. 
     369               z_elem_a( ji,     jpka )       = 0._wp 
     370               z_elem_c( ji,     jpka )       = 0._wp 
     371               z_elem_b( ji,     jpka )       = e3t_abl( jpka ) 
     372               u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk) 
     373            ENDIF 
     374 
     375         END DO 
    328376         !! 
    329377         !! Matrix inversion 
    330378         !! ---------------------------------------------------------- 
    331          DO ji = 2, jpi           
     379         DO ji = 2, jpi 
    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    
     411            DO ji = 1, jpi 
    364412               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               
     413               z_elem_c( ji,     jk ) = -rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal 
    366414               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal 
    367415            END DO 
    368416         END DO 
    369417 
    370          DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)           
     418         DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj) 
    371419            !++ Surface boundary condition 
    372420            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              
     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            IF( ln_topbc_neumann ) THEN 
     434               !++ Top Neumann B.C. 
     435               z_elem_a( ji,     jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
     436               z_elem_c( ji,     jpka ) = 0._wp 
     437               z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka ) 
     438               !v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * v_abl   ( ji, jj, jpka, nt_a ) 
     439            ELSE 
     440               !++ Top Dirichlet B.C. 
     441               z_elem_a( ji,     jpka )       = 0._wp 
     442               z_elem_c( ji,     jpka )       = 0._wp 
     443               z_elem_b( ji,     jpka )       = e3t_abl( jpka ) 
     444               v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk) 
     445            ENDIF 
     446 
     447         END DO 
    394448         !! 
    395449         !! Matrix inversion 
    396450         !! ---------------------------------------------------------- 
    397          DO ji = 1, jpi               
     451         DO ji = 1, jpi 
    398452            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                    
     453            zCF   (ji,   2     ) =   - zcff * z_elem_c( ji,     2       ) 
     454            v_abl (ji,jj,2,nt_a) =     zcff * v_abl   ( ji, jj, 2, nt_a ) 
     455         END DO 
     456 
     457         DO jk = 3, jpka 
     458            DO ji = 1, jpi 
    405459               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) ) 
    406460               zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) 
     
    409463            END DO 
    410464         END DO 
    411              
    412          DO jk = jpkam1,2,-1             
    413             DO ji = 1, jpi     
     465 
     466         DO jk = jpkam1,2,-1 
     467            DO ji = 1, jpi 
    414468               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) 
    415469            END DO 
    416470         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        
     471         ! 
     472      !------------- 
     473      END DO             ! end outer loop 
     474      !------------- 
     475 
     476      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     477      !                            !  5 *** Apply nudging on the dynamics and the tracers 
     478      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     479 
    427480      IF( nn_dyn_restore > 0  ) THEN 
    428          !-------------  
     481         !------------- 
    429482         DO jk = 2, jpka    ! outer loop 
    430          !-------------        
     483         !------------- 
    431484            DO_2D_01_01 
    432485               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) )  
     486               zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) ) 
     487               zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
    435488               zmsk  = msk_abl(ji,jj) 
    436489               zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2   & 
    437490                  &  + jp_alp1_dyn * zsig    + jp_alp0_dyn 
    438491               zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points 
    439                                                              ! rn_Dt = rDt_abl / nn_fsbc                           
     492                                                             ! rn_Dt = rDt_abl / nn_fsbc 
    440493               zcff  = zcff * rest_eq(ji,jj) 
    441                z_cft( ji, jj, jk ) = zcff 
    442494               u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  u_abl( ji, jj, jk, nt_a )   & 
    443                   &                               + zcff   * pu_dta( ji, jj, jk       )                       
     495                  &                               + zcff   * pu_dta( ji, jj, jk       ) 
    444496               v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  v_abl( ji, jj, jk, nt_a )   & 
    445497                  &                               + zcff   * pv_dta( ji, jj, jk       ) 
     
    447499         !------------- 
    448500         END DO             ! end outer loop 
    449          !-------------                
     501         !------------- 
    450502      END IF 
    451503 
    452       !-------------  
     504      !------------- 
    453505      DO jk = 2, jpka    ! outer loop 
    454       !-------------        
     506      !------------- 
    455507         DO_2D_11_11 
    456508            zcff1 = pblh( ji, jj ) 
    457509            zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) ) 
    458             zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )  
     510            zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) ) 
    459511            zmsk  = msk_abl(ji,jj) 
    460512            zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2   & 
    461513               &  + jp_alp1_tra * zsig    + jp_alp0_tra 
    462514            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 
     515                                                          ! rn_Dt = rDt_abl / nn_fsbc 
    465516            tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta )   & 
    466517               &                                       + zcff   * pt_dta( ji, jj, jk ) 
    467              
     518 
    468519            tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa )   & 
    469520               &                                       + zcff   * pq_dta( ji, jj, jk ) 
    470              
     521 
    471522         END_2D 
    472523      !------------- 
    473524      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. ) 
     525      !------------- 
     526      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     527      !                            !  6 *** MPI exchanges 
     528      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     529      ! 
     530      CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a)      , 'T', -1., kfillmode = jpfillnothing ) 
    480531      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... 
    481532      ! 
     533#if defined key_iomput 
    482534      ! first ABL level 
    483535      IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl",   u_abl(:,:,2,nt_a      ) ) 
     
    497549      IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka           ) ) 
    498550      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           ) ) 
     551      IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxlm_abl(:,:,2:jpka          ) ) 
    500552      IF ( iom_use("pblh"   ) ) CALL iom_put (  "pblh"  ,    pblh(:,:                  ) ) 
    501553      ! debug (to be removed) 
     
    504556      IF ( iom_use("t_dta") ) CALL iom_put ( "t_dta",  pt_dta(:,:,2:jpka) ) 
    505557      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) ) 
    507558      IF( ln_geos_winds ) THEN 
    508559         IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2           ) ) 
     
    513564         IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp)   ) 
    514565      END IF 
    515       ! 
     566#endif 
    516567      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    517568      !                            !  7 *** Finalize flux computation 
    518       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  
    519  
     569      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     570      ! 
    520571      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               
     572         ztemp          = tq_abl  ( ji, jj, 2, nt_a, jp_ta ) 
     573         zhumi          = tq_abl  ( ji, jj, 2, nt_a, jp_qa ) 
     574         zcff           = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 
     575         psen( ji, jj ) =    - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )   !GS: negative sign to respect aerobulk convention 
     576         pevp( ji, jj ) = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) ) 
     577         rhoa( ji, jj ) = zcff 
    529578      END_2D 
    530        
     579 
    531580      DO_2D_01_01 
    532          zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji  ,jj) + pssu(ji-1,jj) )   
    533          zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj  ) + pssv(ji,jj-1) )  
     581         zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * rn_vfac * ( pssu(ji,jj) + pssu(ji-1,jj  ) ) 
     582         zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * rn_vfac * ( pssv(ji,jj) + pssv(ji  ,jj-1) ) 
    534583      END_2D 
    535       !  
     584      ! 
    536585      CALL lbc_lnk_multi( 'ablmod', zwnd_i(:,:) , 'T', -1., zwnd_j(:,:) , 'T', -1. ) 
    537586      ! 
    538587      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
    539588      DO_2D_11_11 
    540          zcff          = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   & 
    541             &                 + zwnd_j(ji,jj) * zwnd_j(ji,jj)  )  ! * msk_abl(ji,jj) 
    542          zztmp         = rhoa(ji,jj) * pcd_du(ji,jj) 
    543           
    544          pwndm (ji,jj) =         zcff 
    545          ptaum (ji,jj) = zztmp * zcff 
    546          zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
    547          zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
     589         zcff            = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   & 
     590            &                   + zwnd_j(ji,jj) * zwnd_j(ji,jj)   )  ! * msk_abl(ji,jj) 
     591         zztmp           = rhoa(ji,jj) * pcd_du(ji,jj) 
     592         pwndm (ji,jj)   =         zcff 
     593         ptaum (ji,jj)   = zztmp * zcff 
     594         zwnd_i(ji,jj)   = zztmp * zwnd_i(ji,jj) 
     595         zwnd_j(ji,jj)   = zztmp * zwnd_j(ji,jj) 
     596         taux_abl(ji,jj) = rhoa(ji,jj) * pcd_du(ji,jj) * zwnd_i(ji,jj) 
     597         tauy_abl(ji,jj) = rhoa(ji,jj) * pcd_du(ji,jj) * zwnd_j(ji,jj) 
    548598      END_2D 
    549599      ! ... utau, vtau at U- and V_points, resp. 
     
    574624         ! ------------------------------------------------------------ ! 
    575625         DO_2D_00_00 
    576              
    577             zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 
    578             zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 
    579     
     626 
     627            zztmp1 = 0.5_wp * ( u_abl(ji+1,jj  ,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 
     628            zztmp2 = 0.5_wp * ( v_abl(ji  ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 
     629 
    580630            ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             & 
    581631               &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          & 
     
    592642      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    593643      !                            !  8 *** Swap time indices for the next timestep 
    594       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  
     644      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    595645      nt_n = 1 + MOD( kt  , 2) 
    596646      nt_a = 1 + MOD( kt+1, 2) 
    597       !     
     647      ! 
    598648!--------------------------------------------------------------------------------------------------- 
    599649   END SUBROUTINE abl_stp 
     
    634684      !!                (= Kz dz[Ub] * dz[Un] ) 
    635685      !! --------------------------------------------------------------------- 
    636       INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn  
     686      INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn 
    637687      INTEGER, DIMENSION(1:jpi          )     ::   ikbl 
    638688      REAL(wp)                                ::   zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv 
    639       REAL(wp)                                ::   zdU, zdV, zcff1,zshear,zbuoy,zsig, zustar2 
    640       REAL(wp)                                ::   zdU2,zdV2       
    641       REAL(wp)                                ::   zwndi,zwndj 
     689      REAL(wp)                                ::   zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2 
     690      REAL(wp)                                ::   zdU2, zdV2, zbuoy1, zbuoy2 
     691      REAL(wp)                                ::   zwndi, zwndj 
    642692      REAL(wp), DIMENSION(1:jpi,      1:jpka) ::   zsh2 
    643693      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) ::   zbn2 
    644       REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF        
     694      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF 
    645695      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_a 
    646696      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_b 
     
    648698      LOGICAL                                 ::   ln_Patankar    = .FALSE. 
    649699      LOGICAL                                 ::   ln_dumpvar     = .FALSE. 
    650       LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl  
     700      LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl 
    651701      ! 
    652702      tind  = nt_n 
     
    660710      !------------- 
    661711         ! 
    662          ! Compute vertical shear          
     712         ! Compute vertical shear 
    663713         DO jk = 2, jpkam1 
    664             DO ji = 1,jpi   
    665                zcff        = 1.0_wp / e3w_abl( jk )**2                 
    666                zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2                            
     714            DO ji = 1,jpi 
     715               zcff        = 1.0_wp / e3w_abl( jk )**2 
     716               zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2 
    667717               zdV         = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2 
    668                zsh2(ji,jk) = zdU+zdV 
    669             END DO 
    670          END DO    
     718               zsh2(ji,jk) = zdU+zdV   !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 ) 
     719            END DO 
     720         END DO 
    671721         ! 
    672722         ! Compute brunt-vaisala frequency 
    673723         DO jk = 2, jpkam1 
    674             DO ji = 1,jpi  
    675                zcff  = grav * itvref / e3w_abl( jk )   
     724            DO ji = 1,jpi 
     725               zcff  = grav * itvref / e3w_abl( jk ) 
    676726               zcff1 =  tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk  , tind, jp_ta) 
    677727               zcff2 =  tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa)        & 
     
    679729               zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 )  !<-- zbn2 defined on (2,jpi) 
    680730            END DO 
    681          END DO  
     731         END DO 
    682732         ! 
    683733         ! Terms for the tridiagonal problem 
    684734         DO jk = 2, jpkam1 
    685             DO ji = 1,jpi  
     735            DO ji = 1,jpi 
    686736               zshear       =                 zsh2( ji,     jk )   ! zsh2 is already multiplied by Avm_abl at this point 
    687                zsh2(ji,jk)  = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation           
    688                zbuoy        = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )  
    689                 
     737               zsh2(ji,jk)  = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation 
     738               zbuoy        = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 
     739 
    690740               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 
    691                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           
     741               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 
    692742               IF( (zbuoy + zshear) .gt. 0.) THEN    ! Patankar trick to avoid negative values of TKE 
    693743                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
    694                      &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)     ! diagonal        
    695                   tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )             ! right-hand-side 
     744                     &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk)      ! diagonal 
     745                  tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )               ! right-hand-side 
    696746               ELSE 
    697747                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
    698                      &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)   &  ! diagonal     
    699                      &                     - e3w_abl(jk) * rDt_abl * zbuoy    
    700                   tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl *  zshear )             ! right-hand-side                      
     748                     &                     + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk)   &  ! diagonal 
     749                     &                     - e3w_abl(jk) * rDt_abl * zbuoy 
     750                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl *  zshear )             ! right-hand-side 
    701751               END IF 
    702752            END DO 
    703          END DO   
    704                              
    705          DO ji = 1,jpi    ! vector opt.             
    706             zesrf   =  MAX( 4.63_wp * ustar2(ji,jj), tke_min )            
    707             zetop   = tke_min              
    708             z_elem_a ( ji,     1    ) = 0._wp; z_elem_c ( ji,     1    ) = 0._wp; z_elem_b ( ji,     1    ) = 1._wp                                                    
    709             z_elem_a ( ji,     jpka ) = 0._wp; z_elem_c ( ji,     jpka ) = 0._wp; z_elem_b ( ji,     jpka ) = 1._wp             
    710             tke_abl( ji, jj,    1, nt_a ) = zesrf 
    711             tke_abl( ji, jj, jpka, nt_a ) = zetop  
    712             zbn2(ji,jj,   1) = zbn2( ji,jj,   2)  
    713             zsh2(ji,      1) = zsh2( ji,      2)                                       
    714             zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1)  
     753         END DO 
     754 
     755         DO ji = 1,jpi    ! vector opt. 
     756            zesrf   =  MAX( rn_esfc * ustar2(ji,jj), tke_min ) 
     757            zetop   = tke_min 
     758 
     759            z_elem_a ( ji,        1 ) = 0._wp; z_elem_c ( ji,     1    ) = 0._wp; z_elem_b ( ji,     1    ) = 1._wp 
     760            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 ) 
     761            z_elem_c ( ji,     jpka ) = 0._wp 
     762            z_elem_b ( ji,     jpka ) = e3w_abl(jpka) - z_elem_a(ji, jpka ) 
     763 
     764            tke_abl  ( ji, jj,    1 , nt_a ) = zesrf 
     765            tke_abl  ( ji, jj, jpka , nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n ) 
     766 
     767            zbn2(ji,jj,   1) = zbn2( ji,jj,   2) 
     768            zsh2(ji,      1) = zsh2( ji,      2) 
     769            zbn2(ji,jj,jpka) = zbn2( ji,jj,jpkam1) 
    715770            zsh2(ji,   jpka) = zsh2( ji  , jpkam1) 
    716          END DO         
     771         END DO 
    717772         !! 
    718773         !! Matrix inversion 
     
    720775         DO ji = 1,jpi 
    721776            zcff                  =  1._wp / z_elem_b( ji, 1 ) 
    722             zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       )  
    723             tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a )  
    724          END DO              
    725  
    726          DO jk = 2, jpka             
     777            zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       ) 
     778            tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a ) 
     779         END DO 
     780 
     781         DO jk = 2, jpka 
    727782            DO ji = 1,jpi 
    728783               zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) ) 
     
    732787            END DO 
    733788         END DO 
    734              
    735          DO jk = jpkam1,1,-1             
     789 
     790         DO jk = jpkam1,1,-1 
    736791            DO ji = 1,jpi 
    737792               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) 
    738793            END DO 
    739794         END DO 
    740           
    741 !!FL should not be needed because of Patankar procedure   
     795 
     796!!FL should not be needed because of Patankar procedure 
    742797         tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min ) 
    743798 
     
    745800         !! Diagnose PBL height 
    746801         !! ---------------------------------------------------------- 
    747           
    748           
    749          !                                                        
     802 
     803 
     804         ! 
    750805         ! arrays zRH, zFC and zCF are available at this point 
    751806         ! and zFC(:, 1 ) = 0. 
    752807         ! diagnose PBL height based on zsh2 and zbn2 
    753808         zFC (  :  ,1) = 0._wp 
    754          ikbl( 1:jpi ) = 0  
    755           
     809         ikbl( 1:jpi ) = 0 
     810 
    756811         DO jk = 2,jpka 
    757             DO ji = 1, jpi  
     812            DO ji = 1, jpi 
    758813               zcff  = ghw_abl( jk-1 ) 
    759814               zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) ) 
     
    781836            ELSE 
    782837               pblh( ji, jj ) = ghw_abl(jpka) 
    783             END IF  
    784          END DO 
    785       !-------------       
    786       END DO       
    787       !------------- 
    788       ! 
    789       ! Optional : could add pblh smoothing if pblh is noisy horizontally ...  
     838            END IF 
     839         END DO 
     840      !------------- 
     841      END DO 
     842      !------------- 
     843      ! 
     844      ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 
    790845      IF(ln_smth_pblh) THEN 
    791          CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) 
     846         CALL lbc_lnk( 'ablmod', pblh, 'T', 1., kfillmode = jpfillnothing) 
    792847         CALL smooth_pblh( pblh, msk_abl ) 
    793          CALL lbc_lnk( 'ablmod', pblh, 'T', 1.)    
     848         CALL lbc_lnk( 'ablmod', pblh, 'T', 1., kfillmode = jpfillnothing) 
    794849      ENDIF 
    795850      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    799854      SELECT CASE ( nn_amxl ) 
    800855      ! 
    801       CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom          
    802 #   define zlup zRH       
    803 #   define zldw zFC  
     856      CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom 
     857#   define zlup zRH 
     858#   define zldw zFC 
    804859         DO jj = 1, jpj     ! outer loop 
    805860            ! 
    806861            DO ji = 1, jpi 
    807                mxl_abl ( ji, jj,    1 )  = 0._wp 
    808                mxl_abl ( ji, jj, jpka )  = mxl_min 
    809                zldw( ji,        1 )  = 0._wp               
    810                zlup( ji,     jpka )  = 0._wp 
    811             END DO 
    812             ! 
    813             DO jk = 2, jpkam1    
    814                DO ji = 1, jpi   
    815                   zbuoy             = MAX( zbn2(ji, jj, jk), rsmall )   
    816                   mxl_abl( ji, jj, jk ) = MAX( mxl_min,  & 
    817                      &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy )   )  
    818                END DO 
    819             END DO    
     862               mxld_abl ( ji, jj,    1 ) = 0._wp 
     863               mxld_abl ( ji, jj, jpka ) = mxl_min 
     864               mxlm_abl ( ji, jj,    1 ) = 0._wp 
     865               mxlm_abl ( ji, jj, jpka ) = mxl_min 
     866               zldw     ( ji,        1 ) = zrough(ji,jj) * rn_Lsfc 
     867               zlup     ( ji,     jpka ) = mxl_min 
     868            END DO 
     869            ! 
     870            DO jk = 2, jpkam1 
     871               DO ji = 1, jpi 
     872                  zbuoy             = MAX( zbn2(ji, jj, jk), rsmall ) 
     873                  mxlm_abl( ji, jj, jk ) = MAX( mxl_min,  & 
     874                     &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy )   ) 
     875               END DO 
     876            END DO 
    820877            ! 
    821878            ! Limit mxl 
    822             DO jk = jpkam1,1,-1    
    823                DO ji = 1, jpi  
    824                   zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxl_abl(ji, jj, jk) ) 
    825                END DO 
    826             END DO    
     879            DO jk = jpkam1,1,-1 
     880               DO ji = 1, jpi 
     881                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 
     882               END DO 
     883            END DO 
    827884            ! 
    828885            DO jk = 2, jpka 
    829                DO ji = 1, jpi   
    830                   zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxl_abl(ji, jj, jk) )    
    831                END DO 
    832             END DO  
     886               DO ji = 1, jpi 
     887                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 
     888               END DO 
     889            END DO 
     890            ! 
     891!            DO jk = 1, jpka 
     892!               DO ji = 1, jpi 
     893!                  mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     894!                  mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ),  zlup( ji, jk ) ) 
     895!               END DO 
     896!            END DO 
    833897            ! 
    834898            DO jk = 1, jpka 
    835899               DO ji = 1, jpi 
    836                   mxl_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )    
    837                END DO 
    838             END DO                          
    839             ! 
    840          END DO 
    841 #   undef zlup       
    842 #   undef zldw   
    843          ! 
    844          ! 
    845       CASE ( 1 )             ! length-scale computed as the distance to the PBL height 
    846          DO jj = 1,jpj      ! outer loop 
    847             ! 
    848             DO ji = 1, jpi   ! vector opt. 
    849                zcff = 1._wp / pblh( ji, jj )     ! inverse of hbl               
    850                DO jk = 1, jpka               
    851                   zsig  = MIN( zcff * ghw_abl( jk ), 1. )     
    852                   zcff1 = pblh( ji, jj )                  
    853                   mxl_abl( ji, jj, jk ) =  mxl_min                           & 
    854                      &  +  zsig         * (  amx1*zcff1 + bmx1*mxl_min ) & 
    855                      &  +  zsig *  zsig * (  amx2*zcff1 + bmx2*mxl_min ) & 
    856                      &  +  zsig**3      * (  amx3*zcff1 + bmx3*mxl_min ) & 
    857                      &  +  zsig**4      * (  amx4*zcff1 + bmx4*mxl_min ) & 
    858                      &  +  zsig**5      * (  amx5*zcff1 + bmx5*mxl_min ) 
    859                END DO 
    860             END DO  
    861             !             
    862          END DO             
    863          ! 
     900                  zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     901!                  zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp) 
     902                  mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 
     903                  !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ),  zlup( ji, jk ) ) 
     904                  mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min ) 
     905               END DO 
     906            END DO 
     907            ! 
     908         END DO 
     909#   undef zlup 
     910#   undef zldw 
     911         ! 
     912         ! 
     913      CASE ( 1 )           ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom 
     914#   define zlup zRH 
     915#   define zldw zFC 
     916         DO jj = 1, jpj     ! outer loop 
     917            ! 
     918            DO jk = 2, jpkam1 
     919               DO ji = 1,jpi 
     920                              zcff        = 1.0_wp / e3w_abl( jk )**2 
     921                  zdU         = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2 
     922                  zdV         = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2 
     923                  zsh2(ji,jk) = SQRT(zdU+zdV)   !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 
     924                           ENDDO 
     925                        ENDDO 
     926                        ! 
     927            DO ji = 1, jpi 
     928               zcff                      = zrough(ji,jj) * rn_Lsfc 
     929               mxld_abl ( ji, jj,    1 ) = zcff 
     930               mxld_abl ( ji, jj, jpka ) = mxl_min 
     931               mxlm_abl ( ji, jj,    1 ) = zcff 
     932               mxlm_abl ( ji, jj, jpka ) = mxl_min 
     933               zldw     ( ji,        1 ) = zcff 
     934               zlup     ( ji,     jpka ) = mxl_min 
     935            END DO 
     936            ! 
     937            DO jk = 2, jpkam1 
     938               DO ji = 1, jpi 
     939                  zbuoy    = MAX( zbn2(ji, jj, jk), rsmall ) 
     940                  zcff     = 2.*SQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rod*zsh2(ji,jk) & 
     941                                &             + SQRT( rn_Rod*rn_Rod*zsh2(ji,jk)*zsh2(ji,jk)+2.*zbuoy ) ) 
     942                                  mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff ) 
     943               END DO 
     944            END DO 
     945            ! 
     946            ! Limit mxl 
     947            DO jk = jpkam1,1,-1 
     948               DO ji = 1, jpi 
     949                  zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) ) 
     950               END DO 
     951            END DO 
     952            ! 
     953            DO jk = 2, jpka 
     954               DO ji = 1, jpi 
     955                  zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) ) 
     956               END DO 
     957            END DO 
     958            ! 
     959            DO jk = 1, jpka 
     960               DO ji = 1, jpi 
     961                  !mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     962                  !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp) 
     963                  zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     964                  mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 
     965                  !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) ) 
     966                  mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min ) 
     967               END DO 
     968            END DO 
     969 ! 
     970         END DO 
     971#   undef zlup 
     972#   undef zldw 
     973!         ! 
    864974      CASE ( 2 )           ! Bougeault & Lacarrere 89 length-scale 
    865975         ! 
    866 #   define zlup zRH       
    867 #   define zldw zFC            
     976#   define zlup zRH 
     977#   define zldw zFC 
    868978! zCF is used for matrix inversion 
    869 !              
     979! 
    870980       DO jj = 1, jpj      ! outer loop 
    871           
    872          DO ji = 1, jpi            
    873             zlup( ji,    1 ) = mxl_min 
    874             zldw( ji,    1 ) = mxl_min                
     981 
     982         DO ji = 1, jpi 
     983            zcff             = zrough(ji,jj) * rn_Lsfc 
     984            zlup( ji,    1 ) = zcff 
     985            zldw( ji,    1 ) = zcff 
    875986            zlup( ji, jpka ) = mxl_min 
    876             zldw( ji, jpka ) = mxl_min                 
    877          END DO             
    878           
     987            zldw( ji, jpka ) = mxl_min 
     988         END DO 
     989 
    879990         DO jk = 2,jpka-1 
    880991            DO ji = 1, jpi 
    881992               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk) 
    882                zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)         
    883             END DO 
    884          END DO          
     993               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1) 
     994            END DO 
     995         END DO 
    885996         !! 
    886997         !! BL89 search for lup 
    887          !! ----------------------------------------------------------            
    888          DO jk=2,jpka-1  
     998         !! ---------------------------------------------------------- 
     999         DO jk=2,jpka-1 
    8891000            ! 
    8901001            DO ji = 1, jpi 
     
    8921003               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a ) 
    8931004               ln_foundl(ji ) = .false. 
    894             END DO    
    895             !            
     1005            END DO 
     1006            ! 
    8961007            DO jkup=jk+1,jpka-1 
    8971008               DO ji = 1, jpi 
     1009                  zbuoy1             = MAX( zbn2(ji,jj,jkup  ), rsmall ) 
     1010                  zbuoy2             = MAX( zbn2(ji,jj,jkup-1), rsmall ) 
    8981011                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * & 
    899                      &               ( zbn2(ji,jj,jkup  )*(ghw_abl(jkup  )-ghw_abl(jk)) & 
    900                      &               + zbn2(ji,jj,jkup-1)*(ghw_abl(jkup-1)-ghw_abl(jk)) ) 
     1012                     &               ( zbuoy1*(ghw_abl(jkup  )-ghw_abl(jk)) & 
     1013                     &               + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) ) 
    9011014                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 
    9021015                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk) 
    903                      zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)                    
     1016                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 
    9041017                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   & 
    905                         &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) )  
    906                      zlup(ji,jk) = zcff                 
     1018                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) ) 
     1019                     zlup(ji,jk) = zcff 
     1020                     zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk) 
    9071021                     ln_foundl(ji) = .true. 
    9081022                  END IF 
     
    9101024            END DO 
    9111025            ! 
    912          END DO    
     1026         END DO 
    9131027         !! 
    9141028         !! BL89 search for ldwn 
    915          !! ----------------------------------------------------------           
    916          DO jk=2,jpka-1          
     1029         !! ---------------------------------------------------------- 
     1030         DO jk=2,jpka-1 
    9171031            ! 
    9181032            DO ji = 1, jpi 
     
    9201034               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a ) 
    9211035               ln_foundl(ji ) = .false. 
    922             END DO   
    923             !    
     1036            END DO 
     1037            ! 
    9241038            DO jkdwn=jk-1,1,-1 
    925                DO ji = 1, jpi              
     1039               DO ji = 1, jpi 
     1040                  zbuoy1             = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 
     1041                  zbuoy2             = MAX( zbn2(ji,jj,jkdwn  ), rsmall ) 
    9261042                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  & 
    927                      &               * ( zbn2(ji,jj,jkdwn+1)*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 
    928                                        + zbn2(ji,jj,jkdwn  )*(ghw_abl(jk)-ghw_abl(jkdwn  )) )    
    929                   IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN  
     1043                     &               * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 
     1044                                       + zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn  )) ) 
     1045                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN 
    9301046                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 
    931                      zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )               
     1047                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  ) 
    9321048                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   & 
    933                         &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) )  
    934                      zldw(ji,jk) = zcff  
    935                      ln_foundl(ji) = .true.                
    936                   END IF                                    
    937                END DO 
    938             END DO 
    939             !       
     1049                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) ) 
     1050                     zldw(ji,jk) = zcff 
     1051                     zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  ) 
     1052                     ln_foundl(ji) = .true. 
     1053                  END IF 
     1054               END DO 
     1055            END DO 
     1056            ! 
    9401057         END DO 
    9411058 
    9421059         DO jk = 1, jpka 
    943             DO ji = 1, jpi   
    944                mxl_abl( ji, jj, jk ) = MAX( SQRT( zldw( ji, jk ) * zlup( ji, jk ) ), mxl_min )   
    945             END DO 
    946          END DO   
     1060            DO ji = 1, jpi 
     1061               !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp) 
     1062               zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     1063               mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 
     1064               mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min ) 
     1065            END DO 
     1066         END DO 
    9471067 
    9481068      END DO 
    949 #   undef zlup       
    950 #   undef zldw            
    951          ! 
    952       END SELECT       
     1069#   undef zlup 
     1070#   undef zldw 
     1071         ! 
     1072     CASE ( 3 )           ! Bougeault & Lacarrere 89 length-scale 
     1073         ! 
     1074#   define zlup zRH 
     1075#   define zldw zFC 
     1076! zCF is used for matrix inversion 
     1077! 
     1078       DO jj = 1, jpj      ! outer loop 
     1079          ! 
     1080          DO jk = 2, jpkam1 
     1081             DO ji = 1,jpi 
     1082                            zcff        = 1.0_wp / e3w_abl( jk )**2 
     1083                zdU         = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2 
     1084                zdV         = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2 
     1085                zsh2(ji,jk) = SQRT(zdU+zdV)   !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 ) 
     1086                         ENDDO 
     1087                  ENDDO 
     1088          zsh2(:,      1) = zsh2( :,      2) 
     1089          zsh2(:,   jpka) = zsh2( :, jpkam1) 
     1090 
     1091                 DO ji = 1, jpi 
     1092            zcff              = zrough(ji,jj) * rn_Lsfc 
     1093                        zlup( ji,    1 )  = zcff 
     1094            zldw( ji,    1 )  = zcff 
     1095            zlup( ji, jpka ) = mxl_min 
     1096            zldw( ji, jpka ) = mxl_min 
     1097         END DO 
     1098 
     1099         DO jk = 2,jpka-1 
     1100            DO ji = 1, jpi 
     1101               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk) 
     1102               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1) 
     1103            END DO 
     1104         END DO 
     1105         !! 
     1106         !! BL89 search for lup 
     1107         !! ---------------------------------------------------------- 
     1108         DO jk=2,jpka-1 
     1109            ! 
     1110            DO ji = 1, jpi 
     1111               zCF(ji,1:jpka) = 0._wp 
     1112               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a ) 
     1113               ln_foundl(ji ) = .false. 
     1114            END DO 
     1115            ! 
     1116            DO jkup=jk+1,jpka-1 
     1117               DO ji = 1, jpi 
     1118                  zbuoy1             = MAX( zbn2(ji,jj,jkup  ), rsmall ) 
     1119                  zbuoy2             = MAX( zbn2(ji,jj,jkup-1), rsmall ) 
     1120                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) *                 & 
     1121                     &               ( zbuoy1*(ghw_abl(jkup  )-ghw_abl(jk))      & 
     1122                     &               + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) )    & 
     1123                                         &                            + 0.5_wp * e3t_abl(jkup) * rn_Rod *        & 
     1124                                         &               ( SQRT(tke_abl( ji, jj, jkup  , nt_a ))*zsh2(ji,jkup  ) & 
     1125                                         &               + SQRT(tke_abl( ji, jj, jkup-1, nt_a ))*zsh2(ji,jkup-1) ) 
     1126 
     1127                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN 
     1128                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk) 
     1129                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk) 
     1130                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   & 
     1131                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) ) 
     1132                     zlup(ji,jk) = zcff 
     1133                     zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk) 
     1134                     ln_foundl(ji) = .true. 
     1135                  END IF 
     1136               END DO 
     1137            END DO 
     1138            ! 
     1139         END DO 
     1140         !! 
     1141         !! BL89 search for ldwn 
     1142         !! ---------------------------------------------------------- 
     1143         DO jk=2,jpka-1 
     1144            ! 
     1145            DO ji = 1, jpi 
     1146               zCF(ji,1:jpka) = 0._wp 
     1147               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a ) 
     1148               ln_foundl(ji ) = .false. 
     1149            END DO 
     1150            ! 
     1151            DO jkdwn=jk-1,1,-1 
     1152               DO ji = 1, jpi 
     1153                  zbuoy1             = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 
     1154                  zbuoy2             = MAX( zbn2(ji,jj,jkdwn  ), rsmall ) 
     1155                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  & 
     1156                     &               * (zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1))    & 
     1157                     &                 +zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn  )) )  & 
     1158                                         &                            + 0.5_wp * e3t_abl(jkup) * rn_Rod *          & 
     1159                                         &               ( SQRT(tke_abl( ji, jj, jkdwn+1, nt_a ))*zsh2(ji,jkdwn+1) & 
     1160                                         &               + SQRT(tke_abl( ji, jj, jkdwn  , nt_a ))*zsh2(ji,jkdwn  ) ) 
     1161 
     1162                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN 
     1163                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1) 
     1164                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  ) 
     1165                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   & 
     1166                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) ) 
     1167                     zldw(ji,jk) = zcff 
     1168                                         zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  ) 
     1169                     ln_foundl(ji) = .true. 
     1170                  END IF 
     1171               END DO 
     1172            END DO 
     1173            ! 
     1174         END DO 
     1175 
     1176         DO jk = 1, jpka 
     1177            DO ji = 1, jpi 
     1178               !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp) 
     1179               zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) ) 
     1180               mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 
     1181               mxld_abl( ji, jj, jk ) = MAX(  MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min ) 
     1182            END DO 
     1183         END DO 
     1184 
     1185      END DO 
     1186#   undef zlup 
     1187#   undef zldw 
     1188         ! 
     1189         ! 
     1190         ! 
     1191      END SELECT 
    9531192      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    9541193      !                            !  Finalize the computation of turbulent visc./diff. 
    9551194      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    956        
     1195 
    9571196      !------------- 
    9581197      DO jj = 1, jpj     ! outer loop 
    9591198      !------------- 
    960          DO jk = 1, jpka    
     1199         DO jk = 1, jpka 
    9611200            DO ji = 1, jpi  ! vector opt. 
    962                zcff              = MAX( rn_phimax, rn_Ric * mxl_abl( ji, jj, jk ) * mxl_abl( ji, jj, jk )  & 
    963                   &                                       * zbn2(ji, jj, jk) / tke_abl( ji, jj, jk, nt_a ) )  
    964                zcff2             =  1. / ( 1. + zcff )   !<-- phi_z(z) 
    965                zcff              = mxl_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) ) 
    966 !!FL: MAX function probably useless because of the definition of mxl_min              
     1201               zcff  = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk )  & 
     1202               &     * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) ) 
     1203               zcff2 =  1. / ( 1. + zcff )   !<-- phi_z(z) 
     1204               zcff  = mxlm_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) ) 
     1205               !!FL: MAX function probably useless because of the definition of mxl_min 
    9671206               Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff         , avm_bak   ) 
    968                Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   )                               
    969             END DO 
    970          END DO 
    971       !-------------          
    972       END DO       
     1207               Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   ) 
     1208            END DO 
     1209         END DO 
     1210      !------------- 
     1211      END DO 
    9731212      !------------- 
    9741213 
     
    9881227      !! 
    9891228      !! --------------------------------------------------------------------- 
    990      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: msk    
    991      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d 
     1229      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: msk 
     1230      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d 
    9921231      INTEGER                                     :: ji,jj 
    993      REAL(wp)                                    :: smth_a, smth_b 
    994      REAL(wp), DIMENSION(jpi,jpj)                :: zdX,zdY,zFX,zFY 
    995      REAL(wp)                                    :: zumsk,zvmsk 
     1232      REAL(wp)                                    :: smth_a, smth_b 
     1233      REAL(wp), DIMENSION(jpi,jpj)                :: zdX,zdY,zFX,zFY 
     1234      REAL(wp)                                    :: zumsk,zvmsk 
    9961235      !! 
    9971236      !!========================================================= 
     
    10051244         zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji  ,jj ) ) * zumsk 
    10061245      END_2D 
    1007        
    1008      DO_2D_10_11 
     1246 
     1247      DO_2D_10_11 
    10091248         zvmsk = msk(ji,jj) * msk(ji,jj+1) 
    10101249         zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji  ,jj ) ) * zvmsk 
    1011      END_2D 
    1012        
    1013      DO_2D_10_00 
     1250      END_2D 
     1251 
     1252      DO_2D_10_00 
    10141253         zFY ( ji, jj  ) =   zdY ( ji, jj   )                        & 
    10151254            & +  smth_a*  ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 ))   & 
    10161255            &            -  (zdX ( ji, jj   ) - zdX( ji-1, jj   ))  ) 
    1017      END_2D 
     1256      END_2D 
    10181257 
    10191258      DO_2D_00_10 
     
    10291268  &                 +zFY( ji, jj ) - zFY( ji, jj-1 )  ) 
    10301269      END_2D 
    1031      !! 
     1270 
    10321271!--------------------------------------------------------------------------------------------------- 
    10331272   END SUBROUTINE smooth_pblh 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablrst.F90

    r11945 r12588  
    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(:,:,:           ) ) 
    120119      CALL iom_rstput( iter, nitrst, numraw,    'pblh',    pblh(:,:             ) ) 
    121120      ! 
     
    172171      CALL iom_get( numrar, jpdom_autoglo, 'avm_abl', avm_abl(:,:,:           ) ) 
    173172      CALL iom_get( numrar, jpdom_autoglo, 'avt_abl', avt_abl(:,:,:           ) ) 
    174       CALL iom_get( numrar, jpdom_autoglo, 'mxl_abl', mxl_abl(:,:,:           ) ) 
    175173      CALL iom_get( numrar, jpdom_autoglo,    'pblh',    pblh(:,:             ) ) 
    176174      CALL iom_delay_rst( 'READ', 'ABL', numrar )   ! read only abl delayed global communication variables 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/par_abl.F90

    r12489 r12588  
    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.  !: smoothing of atmospheric PBL height     
    3031 
    31    CHARACTER(len=256), PUBLIC ::   cn_ablrst_in     !: suffix of abl restart name (input) 
    32    CHARACTER(len=256), PUBLIC ::   cn_ablrst_out    !: suffix of abl restart name (output) 
    33    CHARACTER(len=256), PUBLIC ::   cn_ablrst_indir  !: abl restart input directory 
    34    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 
    3537 
    3638   !!--------------------------------------------------------------------- 
     
    4547   REAL(wp), PUBLIC, PARAMETER ::   rn_Cek    = 258._wp                   !: Ekman constant for Richardson number  
    4648   REAL(wp), PUBLIC, PARAMETER ::   rn_epssfc = 1._wp / ( 1._wp + 2.8_wp * 2.8_wp ) 
    47    REAL(wp), PUBLIC            ::   rn_ceps                       !: namelist parameter 
    48    REAL(wp), PUBLIC            ::   rn_cm                         !: namelist parameter 
    49    REAL(wp), PUBLIC            ::   rn_ct                         !: namelist parameter 
    50    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  
    5153   REAL(wp), PUBLIC            ::   rn_Rod                        !: namelist parameter    
    5254   REAL(wp), PUBLIC            ::   rn_Sch     
     55   REAL(wp), PUBLIC            ::   rn_Esfc 
     56   REAL(wp), PUBLIC            ::   rn_Lsfc 
    5357   REAL(wp), PUBLIC            ::   mxl_min     
    5458   REAL(wp), PUBLIC            ::   rn_ldyn_min                   !: namelist parameter 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/sbcabl.F90

    r12549 r12588  
    6868      LOGICAL            ::   lluldl 
    6969      NAMELIST/namsbc_abl/ cn_dir, cn_dom, cn_ablrst_in, cn_ablrst_out,           & 
    70          &                 cn_ablrst_indir, cn_ablrst_outdir,                     & 
     70         &                 cn_ablrst_indir, cn_ablrst_outdir, ln_rstart_abl,      & 
    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 ) 
    167  
     168      rn_Esfc =  1._wp / SQRT(rn_cm*rn_ceps) 
     169      rn_Lsfc = vkarmn * SQRT(SQRT(rn_cm*rn_ceps)) / rn_cm 
     170      
    168171      IF(lwp) THEN 
    169172         WRITE(numout,*) 
     
    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==1) 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 
     
    248255         zcff         = 2._wp * omega * SIN( rad * 90._wp )   !++ fmax 
    249256         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 
    253257      ELSE 
    254258         rest_eq(:,:) = 1._wp 
     
    267271 
    268272      ! initialize ABL from data or restart 
    269       IF( ln_rstart ) THEN 
     273      IF( ln_rstart_abl ) THEN 
    270274         CALL abl_rst_read 
    271275      ELSE 
    272276         CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step 
    273277 
    274          u_abl(:,:,:,nt_n      ) = sf(jp_wndi)%fnow(:,:,:) 
    275          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(:,:,:) 
    276280         tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:) 
    277281         tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:) 
     
    280284         avm_abl(:,:,:          ) = avm_bak 
    281285         avt_abl(:,:,:          ) = avt_bak 
    282          mxl_abl(:,:,:          ) = mxl_min 
    283286         pblh   (:,:            ) = ghw_abl( 3 )  !<-- assume that the pbl contains 3 grid points 
    284287         u_abl  (:,:,:,nt_a     ) = 0._wp 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ICE/icesbc.F90

    r12377 r12588  
    148148      CASE( jp_blk, jp_abl )  !--- bulk formulation & ABL formulation 
    149149                                  CALL blk_ice_2    ( t_su, h_s, h_i, alb_ice, sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),    & 
    150             &                                           sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) )    !  
     150            &                                           sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) ) 
    151151         IF( ln_mixcpl        )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i ) 
    152152         IF( nn_flxdist /= -1 )   CALL ice_flx_dist   ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist ) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90

    r12489 r12588  
    104104   REAL(wp)         ::   rn_zu     ! z(u)   : height of wind measurements 
    105105   ! 
    106    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice 
    107    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme) 
    108    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 
     106   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme and ABL) 
     107   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice 
     108   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 
    109109 
    110110   LOGICAL  ::   ln_skin_cs     ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB 
     
    113113   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 
    114114   LOGICAL  ::   ln_humi_rlh    ! humidity read in files ("sn_humi") is relative humidity     [%] if .true. !LB 
     115   LOGICAL  ::   ln_tpot        !!GS: flag to compute or not potential temperature 
    115116   ! 
    116117   INTEGER  ::   nhumi          ! choice of the bulk algorithm 
     
    170171         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm 
    171172         &                 cn_dir , rn_zqt, rn_zu,                                    & 
    172          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           & 
     173         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, ln_tpot,  & 
    173174         &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh  ! cool-skin / warm-layer !LB 
    174175      !!--------------------------------------------------------------------- 
     
    593594         !#LB: because AGRIF hates functions that return something else than a scalar, need to 
    594595         !     use scalar version of gamma_moist() ... 
    595          DO_2D_11_11 
    596             ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 
    597          END_2D 
     596         IF( ln_tpot ) THEN 
     597            DO_2D_11_11 
     598               ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 
     599            END_2D 
     600         ELSE 
     601            ztpot = ptair(:,:) 
     602         ENDIF 
    598603      ENDIF 
    599604 
     
    663668      ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation 
    664669         CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 
    665             &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         & 
    666             &               wndm(:,:), zU_zu(:,:), pslp(:,:),                 & 
    667             &               taum(:,:), psen(:,:), zqla(:,:),                  & 
    668             &               pEvap=pevp(:,:), prhoa=rhoa(:,:) ) 
     670            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),          & 
     671            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                  & 
     672            &               taum(:,:), psen(:,:), zqla(:,:),                   & 
     673            &               pEvap=pevp(:,:), prhoa=rhoa(:,:)                   ) 
    669674 
    670675         zqla(:,:) = zqla(:,:) * tmask(:,:,1) 
     
    893898 
    894899      ! local scalars ( place there for vector optimisation purposes) 
    895       !IF (ln_abl) rhoa  (:,:)  = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI) 
    896900      zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 
    897901 
     
    913917         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   & 
    914918            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' ) 
    915       ELSE 
     919      ELSE ! ln_abl 
    916920         zztmp1 = 11637800.0_wp 
    917921         zztmp2 =    -5897.8_wp 
Note: See TracChangeset for help on using the changeset viewer.