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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 12749 for NEMO – NEMO

Changeset 12749 for NEMO


Ignore:
Timestamp:
2020-04-15T15:49:42+02:00 (4 years ago)
Author:
gsamson
Message:

ABL code cleaning and small bugfixes; modify shared and ORCA2_ICE_ABL cfgs files accordingly (ticket #2419)

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

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-oce.xml

    r12063 r12749  
    5656     <field field_ref="t_abl" /> 
    5757     <field field_ref="q_abl" /> 
     58     <field field_ref="uvz1_abl" /> 
     59     <field field_ref="tz1_abl" /> 
     60     <field field_ref="qz1_abl" /> 
     61     <field field_ref="uvz1_dta" /> 
     62     <field field_ref="tz1_dta" /> 
     63     <field field_ref="qz1_dta" /> 
    5864     <field field_ref="pblh" /> 
    5965   </file> 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/cfgs/ORCA2_ICE_ABL/EXPREF/namelist_cfg

    r12489 r12749  
    2222   cn_exp      =  "ORCA2"  !  experience name 
    2323   nn_it000    =       1   !  first time step 
    24    nn_itend    =      16   !  last  time step (std 5475) 
     24   nn_itend    =    5840   !  last  time step (std 5840) 
    2525   nn_date0    =  20130101 !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1) 
    2626   nn_write    =       4   !  frequency of write in the output file   (modulo referenced to nn_it000) 
     
    3131&namdom        !   time and space domain 
    3232!----------------------------------------------------------------------- 
    33    rn_Dt      = 5400.     !  time step for the dynamics and tracer 
     33   rn_Dt       = 5400.     !  time step for the dynamics and tracer 
    3434/ 
    3535!----------------------------------------------------------------------- 
     
    5656   sn_sal = 'data_1m_salinity_nomask'             ,     -1.     ,'vosaline',   .true.    , .true. , 'yearly'  ,    ''    ,    ''    ,    '' 
    5757/ 
     58 
    5859!!====================================================================== 
    5960!!            ***  Surface Boundary Condition namelists  ***          !! 
     
    6970!!   namsbc_rnf      river runoffs                                      (ln_rnf     =T) 
    7071!!   namsbc_apr      Atmospheric Pressure                               (ln_apr_dyn =T) 
    71 !!   namsbc_isf      ice shelf melting/freezing                         (ln_isfcav  =T : read (ln_read_cfg=T) or set or usr_def_zgr ) 
    72 !!   namsbc_iscpl    coupling option between land ice model and ocean   (ln_isfcav  =T) 
    7372!!   namsbc_wave     external fields from wave model                    (ln_wave    =T) 
    7473!!   namberg         iceberg floats                                     (ln_icebergs=T) 
     
    7978!----------------------------------------------------------------------- 
    8079   nn_fsbc     = 1         !  frequency of SBC module call 
    81                            !     (also = the frequency of sea-ice & iceberg model call) 
    82                      ! Type of air-sea fluxes  
     80      !                    !  (control sea-ice & iceberg model call) 
     81                     ! Type of air-sea fluxes 
    8382   ln_blk      = .false.   !  Bulk formulation                          (T => fill namsbc_blk ) 
    8483   ln_abl      = .true.    !  ABL  formulation                          (T => fill namsbc_abl ) 
    8584                     ! Sea-ice : 
    86    nn_ice      = 2         !  =2 or 3 automatically for SI3 or CICE    ("key_si3" or "key_cice") 
    87                            !          except in AGRIF zoom where it has to be specified 
    88                      ! Misc. options of sbc :  
     85   nn_ice      = 2         !  =0 no ice boundary condition 
     86      !                    !  =1 use observed ice-cover                 (  => fill namsbc_iif ) 
     87      !                    !  =2 or 3 for SI3 and CICE, respectively 
     88   ln_ice_embd = .false.   !  =T embedded sea-ice (pressure + mass and salt exchanges) 
     89      !                    !  =F levitating ice (no pressure, mass and salt exchanges) 
     90                     ! Misc. options of sbc : 
    8991   ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr) 
     92   ln_dm2dc    = .true.    !  daily mean to diurnal cycle on short wave 
    9093   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr) 
    91    ln_dm2dc    = .true.    !  daily mean to diurnal cycle on short wave 
     94   nn_fwb      = 2         !  FreshWater Budget: =0 unchecked 
     95      !                    !     =1 global mean of e-p-r set to zero at each time step 
     96      !                    !     =2 annual global mean of e-p-r set to zero 
    9297   ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf) 
    93    nn_fwb      = 2         !  FreshWater Budget:  
    94    !                       !    =2 annual global mean of e-p-r set to zero 
    9598   ln_wave     = .false.   !  Activate coupling with wave  (T => fill namsbc_wave) 
    9699   ln_cdgw     = .false.   !  Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave) 
    97    ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave)  
     100   ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave) 
    98101   nn_sdrift   =  0        !  Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift 
    99102      !                    !   = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 
     
    111114   ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    112115   ln_COARE_3p6 = .false.   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
    113    ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 31) 
    114       rn_zqt      = 10.     !  Air temperature & humidity reference height (m) 
    115       rn_zu       = 10.     !  Wind vector reference height (m) 
     116   ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 45r1) 
    116117      ! 
    117       ! Skin is ONLY available in ECMWF and COARE algorithms: 
    118       ln_skin_cs = .false.  !  use the cool-skin parameterization => set nn_fsbc=1 and ln_dm2dc=.true.! 
    119       ln_skin_wl = .false.  !  use the warm-layer        "        => set nn_fsbc=1 and ln_dm2dc=.true.! 
    120       ! 
    121       ln_humi_sph = .true.  !  humidity specified below in "sn_humi" is specific humidity     [kg/kg] if .true. 
    122       ln_humi_dpt = .false. !  humidity specified below in "sn_humi" is dew-point temperature   [K]   if .true. 
    123       ln_humi_rlh = .false. !  humidity specified below in "sn_humi" is relative humidity       [%]   if .true. 
     118      ln_humi_sph = .true.  !  humidity "sn_humi" is specific humidity  [kg/kg] 
     119      ln_tpot     = .false. !!GS: compute potential temperature or not 
    124120   ! 
    125    cn_dir = './'  !  root directory for the bulk data location 
     121   cn_dir      = './'      !  root directory for the bulk data location 
    126122   !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________! 
    127123   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ !          weights filename            ! rotation ! land/sea mask ! 
     
    146142   cn_dir         = './'      !  root directory for the location of the ABL grid file 
    147143   cn_dom         = 'dom_cfg_abl_L25Z10.nc' 
     144 
     145   ln_rstart_abl  = .false. 
    148146   ln_hpgls_frc   = .true. 
    149147   ln_geos_winds  = .false. 
    150148   nn_dyn_restore = 1 
    151    rn_ldyn_min   =  7.5       !  magnitude of the nudging on ABL dynamics at the bottom of the ABL   [hour] 
     149 
     150   rn_ldyn_min   =  4.5       !  magnitude of the nudging on ABL dynamics at the bottom of the ABL   [hour] 
    152151   rn_ldyn_max   =  1.5       !  magnitude of the nudging on ABL dynamics at the top of the ABL   [hour] 
    153    rn_ltra_min   =  7.5       !  magnitude of the nudging on ABL tracers  at the bottom of the ABL   [hour] 
     152   rn_ltra_min   =  4.5       !  magnitude of the nudging on ABL tracers  at the bottom of the ABL   [hour] 
    154153   rn_ltra_max   =  1.5       !  magnitude of the nudging on ABL tracers  at the top of the ABL   [hour] 
    155154/ 
     
    201200&namberg       !   iceberg parameters                                   (default: OFF) 
    202201!----------------------------------------------------------------------- 
     202   ln_icebergs = .true.    ! activate iceberg floats (force =F with "key_agrif") 
     203 
     204   cn_dir = './'  !  root directory for the location of drag coefficient files 
     205   !______!___________!___________________!______________!______________!_________!___________!__________!__________!_______________! 
     206   !      ! file name ! frequency (hours) !   variable   ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
     207   !      !           !  (if <0  months)  !     name     !   (logical)  !  (T/F ) ! 'monthly' ! filename ! pairing  ! filename      ! 
     208   sn_icb =  'calving',       -1.         , 'calving'    ,   .true.     , .true.  , 'yearly'  , ''       , ''       , '' 
    203209/ 
    204210!!====================================================================== 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/cfgs/SHARED/field_def_nemo-oce.xml

    r12377 r12749  
    440440          <field id="t_dta"      long_name="DTA potential temperature"     standard_name="dta_theta"      unit="K"        /> 
    441441          <field id="q_dta"      long_name="DTA specific humidity"         standard_name="dta_qspe"       unit="kg/kg"    /> 
    442           <field id="coeft"      long_name="ABL nudging coefficient"       standard_name="coeft"          unit=""         /> 
     442          <field id="u_geo"      long_name="GEO i-horizontal velocity"     standard_name="geo_x_velocity" unit="m/s"      /> 
     443          <field id="v_geo"      long_name="GEO j-horizontal velocity"     standard_name="geo_y_velocity" unit="m/s"      /> 
    443444          <field id="tke_abl"    long_name="ABL turbulent kinetic energy"  standard_name="abl_tke"        unit="m2/s2"    /> 
    444445          <field id="avm_abl"    long_name="ABL turbulent viscosity"       standard_name="abl_avm"        unit="m2/s"     /> 
    445446          <field id="avt_abl"    long_name="ABL turbulent diffusivity"     standard_name="abl_avt"        unit="m2/s"     /> 
    446           <field id="mxl_abl"    long_name="ABL mixing length"             standard_name="abl_mxl"        unit="m"        /> 
     447          <field id="mxlm_abl"   long_name="ABL master mixing length"      standard_name="abl_mxlm"       unit="m"        /> 
     448          <field id="mxld_abl"   long_name="ABL dissipative mixing length" standard_name="abl_mxld"       unit="m"        /> 
    447449   </field_group> 
    448450 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/cfgs/SHARED/namelist_ref

    r12589 r12749  
    259259!----------------------------------------------------------------------- 
    260260   !                    !  bulk algorithm : 
    261    ln_NCAR      = .true.    ! "NCAR"      algorithm   (Large and Yeager 2008) 
     261   ln_NCAR      = .false.   ! "NCAR"      algorithm   (Large and Yeager 2008) 
    262262   ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    263263   ln_COARE_3p6 = .false.   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
     
    271271      rn_pfac    = 1.       !  multipl. factor for precipitation (total & snow) 
    272272      rn_efac    = 1.       !  multipl. factor for evaporation (0. or 1.) 
    273       rn_vfac    = 0.       !  multipl. factor for ocean & ice velocity 
     273      rn_vfac    = 1.       !  multipl. factor for ocean & ice velocity 
    274274      !                     !  used to calculate the wind stress 
    275275      !                     ! (0. => absolute or 1. => relative winds) 
     
    277277      ln_skin_wl = .false.  !  use the warm-layer parameterization 
    278278      !                     !   ==> only available in ECMWF and COARE algorithms 
    279       ln_humi_sph = .true. !  humidity "sn_humi" is specific humidity  [kg/kg] 
     279      ln_humi_sph = .false. !  humidity "sn_humi" is specific humidity  [kg/kg] 
    280280      ln_humi_dpt = .false. !  humidity "sn_humi" is dew-point temperature [K] 
    281281      ln_humi_rlh = .false. !  humidity "sn_humi" is relative humidity     [%] 
     282      ln_tpot     = .true.  !!GS: compute potential temperature or not 
    282283   ! 
    283284   cn_dir      = './'      !  root directory for the bulk data location 
     
    315316                              !                                               = 1 equatorial restoring 
    316317                              !                                               = 2 global restoring 
    317    rn_ldyn_min   =  4.5       !  magnitude of the nudging on ABL dynamics at the bottom of the ABL   [hour] 
    318    rn_ldyn_max   =  1.5       !  magnitude of the nudging on ABL dynamics at the top of the ABL   [hour] 
    319    rn_ltra_min   =  4.5       !  magnitude of the nudging on ABL tracers  at the bottom of the ABL   [hour] 
    320    rn_ltra_max   =  1.5       !  magnitude of the nudging on ABL tracers  at the top of the ABL   [hour] 
     318   rn_ldyn_min   =  4.5       ! dynamics nudging magnitude inside the ABL [hour] (~3 rn_Dt) 
     319   rn_ldyn_max   =  1.5       ! dynamics nudging magnitude above  the ABL [hour] (~1 rn_Dt) 
     320   rn_ltra_min   =  4.5       ! tracers  nudging magnitude inside the ABL [hour] (~3 rn_Dt) 
     321   rn_ltra_max   =  1.5       ! tracers  nudging magnitude above  the ABL [hour] (~1 rn_Dt) 
    321322   nn_amxl       =  0         ! mixing length: = 0 Deardorff 80 length-scale 
    322323                              !                = 1 length-scale based on the distance to the PBL height 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/cfgs/ref_cfgs.txt

    r12377 r12749  
    77ORCA2_OFF_TRC OCE TOP OFF 
    88ORCA2_SAS_ICE OCE ICE NST SAS 
    9 ORCA2_ICE_PISCES OCE TOP ICE NST 
     9ORCA2_ICE_PISCES OCE TOP ICE NST ABL 
    1010ORCA2_ICE_ABL OCE ICE ABL 
    11 ORCA2_SAS_ICE_ABL OCE SAS ICE ABL 
    12 ORCA2_ICE OCE ICE 
    1311SPITZ12 OCE ICE 
    1412WED025 OCE ICE 
    15 eORCA025_ICE OCE ICE 
    16 eORCA025_ICE_ABL OCE ICE ABL 
    17 eORCA025_SAS_ICE_ABL OCE SAS ICE ABL 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/abl.F90

    r12588 r12749  
    3939   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       ::   msk_abl 
    4040   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   
    4541   ! 
    4642   INTEGER , PUBLIC :: nt_n, nt_a       !: now / after indices (equal 1 or 2) 
     
    6864         &      mxld_abl(1:jpi,1:jpj,1:jpka            ), &          
    6965         &      mxlm_abl(1:jpi,1:jpj,1:jpka            ), &  
    70          &      cft_abl (1:jpi,1:jpj,1:jpka            ), &  
    7166         &      fft_abl (1:jpi,1:jpj                   ), & 
    7267         &      pblh    (1:jpi,1:jpj                   ), &          
    73          &      taux_abl(1:jpi,1:jpj                   ), &    
    74          &      tauy_abl(1:jpi,1:jpj                   ), &         
    7568         &      msk_abl (1:jpi,1:jpj                   ), & 
    7669         &      rest_eq (1:jpi,1:jpj                   ), &          
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/ablmod.F90

    r12592 r12749  
    3838 
    3939!=================================================================================================== 
    40    SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, &                            ! in 
     40   SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq,            &     ! in 
    4141              &            pu_dta, pv_dta, pt_dta, pq_dta,    & 
    4242              &            pslp_dta, pgu_dta, pgv_dta,        & 
    43               &            pcd_du, psen, pevp,                &                            ! in/out 
     43              &            pcd_du, psen, pevp,                &     ! in/out 
    4444              &            pwndm, ptaui, ptauj, ptaum         & 
    4545#if defined key_si3 
     
    112112      REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2 
    113113      REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce 
    114       LOGICAL             ::   SemiImp_Cor = .FALSE. 
     114      LOGICAL             ::   SemiImp_Cor = .TRUE. 
    115115      ! 
    116116      !!--------------------------------------------------------------------- 
     
    135135         ustar2(ji,jj) = zzoce 
    136136#endif 
    137          zrough(ji,jj) = ght_abl(2) * EXP( - vkarmn / SQRT(Cdn_oce(ji,jj)) ) !<-- recover the value of z0 from Cdn_oce 
     137         zrough(ji,jj) = ght_abl(2) * EXP( - vkarmn / SQRT( MAX( Cdn_oce(ji,jj), 1.e-4 ) ) ) !<-- recover the value of z0 from Cdn_oce 
    138138      END_2D 
    139139      ! 
     
    162162         DO ji = 1, jpi   ! vector opt. 
    163163            ! 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  ) 
     164            z_elem_a( ji,    2 ) = 0._wp 
     165            z_elem_c( ji,    2 ) = - rDt_abl * Avt_abl( ji, jj,    2 ) / e3w_abl(    2 ) 
    166166            ! Homogeneous Neumann at the top 
    167167            z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
    168168            z_elem_c( ji, jpka ) = 0._wp 
    169             z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka ) 
     169            z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) 
    170170         END DO 
    171171 
     
    173173 
    174174            DO jk = 3, jpkam1 
    175                DO ji = 2, jpim1 
    176                !DO ji = 1,jpi  !!GS: to be checked 
    177                   tq_abl  ( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl  ( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side 
     175               !DO ji = 2, jpim1 
     176               DO ji = 1,jpi  !!GS: to be checked if needed 
     177                  tq_abl( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side 
    178178               END DO 
    179179            END DO 
     
    187187                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj) 
    188188#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 
     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 
    191191                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra ) 
    192192               END DO 
     
    199199                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj) 
    200200#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 
     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 
    203203                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra ) 
    204204               END DO 
     
    208208            !! ---------------------------------------------------------- 
    209209            DO ji = 1,jpi 
    210                zcff                       =  1._wp / z_elem_b( ji, 2 ) 
    211                zCF   (ji,   2           ) = - zcff * z_elem_c( ji, 2 ) 
    212                tq_abl(ji,jj,2,nt_a,jtra) =    zcff * tq_abl(ji,jj,2,nt_a,jtra) 
     210               zcff                            =  1._wp / z_elem_b( ji, 2 ) 
     211               zCF   ( ji,     2             ) = - zcff * z_elem_c( ji, 2 ) 
     212               tq_abl( ji, jj, 2, nt_a, jtra ) =   zcff * tq_abl( ji, jj, 2, nt_a, jtra ) 
    213213            END DO 
    214214 
    215215            DO jk = 3, jpka 
    216216               DO ji = 1,jpi 
    217                   zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) ) 
     217                  zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF( ji, jk-1 ) ) 
    218218                  zCF(ji,jk) = - zcff * z_elem_c( ji, jk ) 
    219219                  tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk  ,nt_a,jtra)   & 
    220                   &                - z_elem_a(ji, jk) * tq_abl(ji,jj,jk-1,nt_a,jtra) ) 
     220                     &           - z_elem_a(ji, jk) *  tq_abl(ji,jj,jk-1,nt_a,jtra) ) 
    221221               END DO 
    222222            END DO 
     
    240240      IF( SemiImp_Cor ) THEN 
    241241 
    242       !------------- 
    243       DO jk = 2, jpka    ! outer loop 
    244       !------------- 
    245          ! 
    246          ! Advance u_abl & v_abl to time n+1 
    247          DO_2D_11_11 
    248             zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl )  ! (f dt)**2 
    249  
    250             u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(  & 
    251                &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n )    & 
    252                &                 +  rDt_abl * fft_abl(ji, jj) * v_abl ( ji , jj  , jk, nt_n ) )  & 
    253                &                               / (1._wp + gamma_Cor*gamma_Cor*zcff) 
    254  
    255             v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(  & 
    256                &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n )   & 
    257                &                 -  rDt_abl * fft_abl(ji, jj) * u_abl ( ji   , jj, jk, nt_n )  ) & 
    258                &                                / (1._wp + gamma_Cor*gamma_Cor*zcff) 
    259          END_2D 
    260          ! 
    261       !------------- 
    262       END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 
    263       !------------- 
    264       ! 
    265       !IF( ln_geos_winds ) THEN 
    266       !   DO jj = 1, jpj    ! outer loop 
    267       !      DO jk = 1, jpka 
    268       !         DO ji = 1, jpi 
    269       !            u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   & 
    270       !               &                      - rDt_abl * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk) 
    271       !            v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   & 
    272       !               &                      + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk) 
    273       !         END DO 
    274       !      END DO 
    275       !   END DO 
    276       !END IF 
    277  
    278       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 
    320             DO jk = 1, jpka 
    321                DO ji = 1, jpi 
    322                   u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_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) 
     242         !------------- 
     243         DO jk = 2, jpka    ! outer loop 
     244         !------------- 
     245            ! 
     246            ! Advance u_abl & v_abl to time n+1 
     247            DO_2D_11_11 
     248               zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl )  ! (f dt)**2 
     249 
     250               u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(                                          & 
     251                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff) * u_abl( ji, jj, jk, nt_n )    & 
     252                  &                     + rDt_abl * fft_abl(ji, jj) * v_abl( ji, jj, jk, nt_n ) )  & 
     253                  &                               / (1._wp + gamma_Cor*gamma_Cor*zcff) 
     254    
     255               v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(                                         & 
     256                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff) * v_abl( ji, jj, jk, nt_n )    & 
     257                  &                     - rDt_abl * fft_abl(ji, jj) * u_abl( ji, jj, jk, nt_n ) )  & 
     258                  &                               / (1._wp + gamma_Cor*gamma_Cor*zcff) 
     259            END_2D 
     260            ! 
     261         !------------- 
     262         END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 
     263         !------------- 
     264         ! 
     265         IF( ln_geos_winds ) THEN 
     266            DO jj = 1, jpj    ! outer loop 
     267               DO jk = 1, jpka 
     268                  DO ji = 1, jpi 
     269                     u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   & 
     270                        &                      - rDt_abl * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk) 
     271                     v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   & 
     272                        &                      + rDt_abl * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk) 
     273                  END DO 
     274               END DO 
     275            END DO 
     276         END IF 
     277         ! 
     278         IF( ln_hpgls_frc ) THEN 
     279            DO jj = 1, jpj    ! outer loop 
     280               DO jk = 1, jpka 
     281                  DO ji = 1, jpi 
     282                     u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk) 
     283                     v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk) 
     284                  ENDDO 
    324285               ENDDO 
    325286            ENDDO 
    326          ENDDO 
    327       END IF 
    328  
     287         END IF 
     288 
     289      ELSE ! SemiImp_Cor = .FALSE. 
     290 
     291         IF( ln_geos_winds ) THEN 
     292 
     293            !------------- 
     294            DO jk = 2, jpka    ! outer loop 
     295            !------------- 
     296               ! 
     297               IF( MOD( kt, 2 ) == 0 ) then 
     298                  ! Advance u_abl & v_abl to time n+1 
     299                  DO jj = 1, jpj 
     300                     DO ji = 1, jpi 
     301                        zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_n ) - pgv_dta(ji  ,jj  ,jk)  ) 
     302                        u_abl( ji, jj, jk, nt_a ) =                u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff 
     303                        zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_a ) - pgu_dta(ji  ,jj  ,jk)  ) 
     304                        v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff ) 
     305                        u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  u_abl( ji, jj, jk, nt_a ) 
     306                     END DO 
     307                  END DO 
     308               ELSE 
     309                  DO jj = 1, jpj 
     310                     DO ji = 1, jpi 
     311                        zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_n ) - pgu_dta(ji  ,jj  ,jk)  ) 
     312                        v_abl( ji, jj, jk, nt_a ) =                v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff 
     313                        zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_a ) - pgv_dta(ji  ,jj  ,jk)  ) 
     314                        u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff ) 
     315                        v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  v_abl( ji, jj, jk, nt_a ) 
     316                     END DO 
     317                  END DO 
     318               END IF 
     319               ! 
     320            !------------- 
     321            END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj) 
     322            !------------- 
     323 
     324         ENDIF ! ln_geos_winds 
     325 
     326      ENDIF ! SemiImp_Cor 
    329327      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    330328      !                            !  4 *** Advance u,v to time n+1 
     
    338336         DO jk = 3, jpkam1 
    339337            DO ji = 1, jpi 
    340                z_elem_a( ji,     jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal 
    341                z_elem_c( ji,     jk ) = - rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal 
    342                z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                             !       diagonal 
     338               z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal 
     339               z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal 
     340               z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )  !       diagonal 
    343341            END DO 
    344342         END DO 
     
    346344         DO ji = 2, jpi   ! boundary conditions   (Avm_abl and pcd_du must be available at ji=jpi) 
    347345            !++ Surface boundary condition 
    348             z_elem_a( ji,     2    ) = 0._wp 
    349             z_elem_c( ji,     2    ) = - rDt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2  ) 
     346            z_elem_a( ji, 2 ) = 0._wp 
     347            z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 
    350348            ! 
    351349            zztmp1  = pcd_du(ji, jj) 
     
    359357            u_abl( ji, jj,    2, nt_a ) =      u_abl( ji, jj,    2, nt_a ) + rDt_abl * zztmp2 
    360358 
    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 
     359            ! idealized test cases only 
     360            !IF( ln_topbc_neumann ) THEN 
     361            !   !++ Top Neumann B.C. 
     362            !   z_elem_a( ji,     jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
     363            !   z_elem_c( ji,     jpka ) = 0._wp 
     364            !   z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka ) 
     365            !   !u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * u_abl   ( ji, jj, jpka, nt_a ) 
     366            !ELSE 
    368367               !++ Top Dirichlet B.C. 
    369368               z_elem_a( ji,     jpka )       = 0._wp 
     
    371370               z_elem_b( ji,     jpka )       = e3t_abl( jpka ) 
    372371               u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk) 
    373             ENDIF 
     372            !ENDIF 
    374373 
    375374         END DO 
     
    377376         !! Matrix inversion 
    378377         !! ---------------------------------------------------------- 
    379          DO ji = 2, jpi 
     378         !DO ji = 2, jpi 
     379         DO ji = 1, jpi  !!GS: TBI 
    380380            zcff                 =   1._wp / z_elem_b( ji, 2 ) 
    381381            zCF   (ji,   2     ) =  - zcff * z_elem_c( ji, 2 ) 
     
    410410         DO jk = 3, jpkam1 
    411411            DO ji = 1, jpi 
    412                z_elem_a( ji,     jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal 
    413                z_elem_c( ji,     jk ) = -rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal 
    414                z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal 
     412               z_elem_a( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal 
     413               z_elem_c( ji, jk ) = -rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal 
     414               z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal 
    415415            END DO 
    416416         END DO 
     
    418418         DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj) 
    419419            !++ Surface boundary condition 
    420             z_elem_a( ji,     2    ) = 0._wp 
    421             z_elem_c( ji,     2    ) = - rDt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2  ) 
     420            z_elem_a( ji, 2 ) = 0._wp 
     421            z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 
    422422            ! 
    423423            zztmp1 = pcd_du(ji, jj) 
     
    431431            v_abl( ji, jj,    2, nt_a ) =         v_abl( ji, jj, 2, nt_a ) + rDt_abl * zztmp2 
    432432 
    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 
     433            ! idealized test cases only 
     434            !IF( ln_topbc_neumann ) THEN 
     435            !   !++ Top Neumann B.C. 
     436            !   z_elem_a( ji,     jpka ) = - rDt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 
     437            !   z_elem_c( ji,     jpka ) = 0._wp 
     438            !   z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka ) 
     439            !   !v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * v_abl   ( ji, jj, jpka, nt_a ) 
     440            !ELSE 
    440441               !++ Top Dirichlet B.C. 
    441442               z_elem_a( ji,     jpka )       = 0._wp 
     
    443444               z_elem_b( ji,     jpka )       = e3t_abl( jpka ) 
    444445               v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk) 
    445             ENDIF 
     446            !ENDIF 
    446447 
    447448         END DO 
     
    528529      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    529530      ! 
    530       CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a)      , 'T', -1., kfillmode = jpfillnothing ) 
     531      CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a)      , 'T', -1.                            ) 
    531532      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... 
    532533      ! 
    533534#if defined key_iomput 
    534       ! first ABL level 
     535      ! 2D & first ABL level 
     536      IF ( iom_use("pblh"   ) ) CALL iom_put (    "pblh",    pblh(:,:             ) ) 
    535537      IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl",   u_abl(:,:,2,nt_a      ) ) 
    536538      IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl",   v_abl(:,:,2,nt_a      ) ) 
     
    541543      IF ( iom_use("tz1_dta") ) CALL iom_put ( "tz1_dta",  pt_dta(:,:,2           ) ) 
    542544      IF ( iom_use("qz1_dta") ) CALL iom_put ( "qz1_dta",  pq_dta(:,:,2           ) ) 
    543       ! all ABL levels 
    544       IF ( iom_use("u_abl"  ) ) CALL iom_put ( "u_abl"  ,   u_abl(:,:,2:jpka,nt_a      ) ) 
    545       IF ( iom_use("v_abl"  ) ) CALL iom_put ( "v_abl"  ,   v_abl(:,:,2:jpka,nt_a      ) ) 
    546       IF ( iom_use("t_abl"  ) ) CALL iom_put ( "t_abl"  ,  tq_abl(:,:,2:jpka,nt_a,jp_ta) ) 
    547       IF ( iom_use("q_abl"  ) ) CALL iom_put ( "q_abl"  ,  tq_abl(:,:,2:jpka,nt_a,jp_qa) ) 
    548       IF ( iom_use("tke_abl") ) CALL iom_put ( "tke_abl", tke_abl(:,:,2:jpka,nt_a      ) ) 
    549       IF ( iom_use("avm_abl") ) CALL iom_put ( "avm_abl", avm_abl(:,:,2:jpka           ) ) 
    550       IF ( iom_use("avt_abl") ) CALL iom_put ( "avt_abl", avm_abl(:,:,2:jpka           ) ) 
    551       IF ( iom_use("mxl_abl") ) CALL iom_put ( "mxl_abl", mxlm_abl(:,:,2:jpka          ) ) 
    552       IF ( iom_use("pblh"   ) ) CALL iom_put (  "pblh"  ,    pblh(:,:                  ) ) 
    553       ! debug (to be removed) 
     545      ! debug 2D 
     546      IF( ln_geos_winds ) THEN 
     547         IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2) ) 
     548         IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2) ) 
     549      END IF 
     550      IF( ln_hpgls_frc ) THEN 
     551         IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo",  pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) ) 
     552         IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp) ) 
     553      END IF 
     554      ! 3D (all ABL levels) 
     555      IF ( iom_use("u_abl"   ) ) CALL iom_put ( "u_abl"   ,    u_abl(:,:,2:jpka,nt_a      ) ) 
     556      IF ( iom_use("v_abl"   ) ) CALL iom_put ( "v_abl"   ,    v_abl(:,:,2:jpka,nt_a      ) ) 
     557      IF ( iom_use("t_abl"   ) ) CALL iom_put ( "t_abl"   ,   tq_abl(:,:,2:jpka,nt_a,jp_ta) ) 
     558      IF ( iom_use("q_abl"   ) ) CALL iom_put ( "q_abl"   ,   tq_abl(:,:,2:jpka,nt_a,jp_qa) ) 
     559      IF ( iom_use("tke_abl" ) ) CALL iom_put ( "tke_abl" ,  tke_abl(:,:,2:jpka,nt_a      ) ) 
     560      IF ( iom_use("avm_abl" ) ) CALL iom_put ( "avm_abl" ,  avm_abl(:,:,2:jpka           ) ) 
     561      IF ( iom_use("avt_abl" ) ) CALL iom_put ( "avt_abl" ,  avt_abl(:,:,2:jpka           ) ) 
     562      IF ( iom_use("mxlm_abl") ) CALL iom_put ( "mxlm_abl", mxlm_abl(:,:,2:jpka           ) ) 
     563      IF ( iom_use("mxld_abl") ) CALL iom_put ( "mxld_abl", mxld_abl(:,:,2:jpka           ) ) 
     564      ! debug 3D 
    554565      IF ( iom_use("u_dta") ) CALL iom_put ( "u_dta",  pu_dta(:,:,2:jpka) ) 
    555566      IF ( iom_use("v_dta") ) CALL iom_put ( "v_dta",  pv_dta(:,:,2:jpka) ) 
     
    557568      IF ( iom_use("q_dta") ) CALL iom_put ( "q_dta",  pq_dta(:,:,2:jpka) ) 
    558569      IF( ln_geos_winds ) THEN 
    559          IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo", pgu_dta(:,:,2           ) ) 
    560          IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", pgv_dta(:,:,2           ) ) 
     570         IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo", pgu_dta(:,:,2:jpka) ) 
     571         IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", pgv_dta(:,:,2:jpka) ) 
    561572      END IF 
    562573      IF( ln_hpgls_frc ) THEN 
    563          IF ( iom_use("uz1_geo") ) CALL iom_put ( "uz1_geo",  pgu_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp)  ) 
    564          IF ( iom_use("vz1_geo") ) CALL iom_put ( "vz1_geo", -pgv_dta(:,:,2)/MAX(fft_abl(:,:),2.5e-5_wp)  ) 
     574         IF ( iom_use("u_geo") ) CALL iom_put ( "u_geo",  pgu_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) ) 
     575         IF ( iom_use("v_geo") ) CALL iom_put ( "v_geo", -pgv_dta(:,:,2:jpka)/MAX( RESHAPE( fft_abl(:,:), (/jpi,jpj,jpka-1/), fft_abl(:,:)), 2.5e-5_wp) ) 
    565576      END IF 
    566577#endif 
     
    570581      ! 
    571582      DO_2D_11_11 
    572          ztemp          = tq_abl  ( ji, jj, 2, nt_a, jp_ta ) 
    573          zhumi          = tq_abl  ( ji, jj, 2, nt_a, jp_qa ) 
     583         ztemp          =  tq_abl( ji, jj, 2, nt_a, jp_ta ) 
     584         zhumi          =  tq_abl( ji, jj, 2, nt_a, jp_qa ) 
    574585         zcff           = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 
    575586         psen( ji, jj ) =    - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )   !GS: negative sign to respect aerobulk convention 
     
    587598      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
    588599      DO_2D_11_11 
    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) 
     600         zcff          = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   & 
     601            &                 + zwnd_j(ji,jj) * zwnd_j(ji,jj) )   ! * msk_abl(ji,jj) 
     602         zztmp         = rhoa(ji,jj) * pcd_du(ji,jj) 
     603 
     604         pwndm (ji,jj) =         zcff 
     605         ptaum (ji,jj) = zztmp * zcff 
     606         zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
     607         zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
    598608      END_2D 
    599609      ! ... utau, vtau at U- and V_points, resp. 
     
    620630 
    621631#if defined key_si3 
    622          ! ------------------------------------------------------------ ! 
    623          !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
    624          ! ------------------------------------------------------------ ! 
    625          DO_2D_00_00 
    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  
    630             ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             & 
    631                &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          & 
    632                &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 
    633             ptauj_ice(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             & 
    634                &                      +    rhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          & 
    635                &         * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 
    636          END_2D 
    637          CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. ) 
    638          ! 
    639          IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   & 
    640             &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' ) 
     632      ! ------------------------------------------------------------ ! 
     633      !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     634      ! ------------------------------------------------------------ ! 
     635      DO_2D_00_00 
     636 
     637         zztmp1 = 0.5_wp * ( u_abl(ji+1,jj  ,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 
     638         zztmp2 = 0.5_wp * ( v_abl(ji  ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 
     639 
     640         ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             & 
     641            &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          & 
     642            &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 
     643         ptauj_ice(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             & 
     644            &                      +    rhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          & 
     645            &         * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 
     646      END_2D 
     647      CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. ) 
     648      ! 
     649      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   & 
     650         &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' ) 
    641651#endif 
    642652      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    688698      REAL(wp)                                ::   zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv 
    689699      REAL(wp)                                ::   zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2 
    690       REAL(wp)                                ::   zdU2, zdV2, zbuoy1, zbuoy2 
     700      REAL(wp)                                ::   zdU2, zdV2, zbuoy1, zbuoy2    ! zbuoy for BL89 
    691701      REAL(wp)                                ::   zwndi, zwndj 
    692702      REAL(wp), DIMENSION(1:jpi,      1:jpka) ::   zsh2 
     
    733743         ! Terms for the tridiagonal problem 
    734744         DO jk = 2, jpkam1 
    735             DO ji = 1,jpi 
    736                zshear       =                 zsh2( ji,     jk )   ! zsh2 is already multiplied by Avm_abl at this point 
    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  
    740                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 
    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 
     745            DO ji = 1, jpi 
     746               zshear      = zsh2( ji, jk )                           ! zsh2 is already multiplied by Avm_abl at this point 
     747               zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation 
     748               zbuoy       = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 
     749 
     750               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 
     751               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 
    742752               IF( (zbuoy + zshear) .gt. 0.) THEN    ! Patankar trick to avoid negative values of TKE 
    743                   z_elem_b( ji,     jk )  = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
    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 
     753                  z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
     754                     &               + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk)    ! diagonal 
     755                  tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl * ( zbuoy + zshear ) )       ! right-hand-side 
    746756               ELSE 
    747                   z_elem_b( ji,     jk )  = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
    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 
     757                  z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
     758                     &               + e3w_abl(jk) * rDt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxld_abl(ji,jj,jk)   &  ! diagonal 
     759                     &               - e3w_abl(jk) * rDt_abl * zbuoy 
     760                  tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl *  zshear )                    ! right-hand-side 
    751761               END IF 
    752762            END DO 
     
    754764 
    755765         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) 
    770             zsh2(ji,   jpka) = zsh2( ji  , jpkam1) 
     766            zesrf = MAX( rn_Esfc * ustar2(ji,jj), tke_min ) 
     767            zetop = tke_min 
     768 
     769            z_elem_a ( ji,     1       ) = 0._wp 
     770            z_elem_c ( ji,     1       ) = 0._wp 
     771            z_elem_b ( ji,     1       ) = 1._wp 
     772            tke_abl  ( ji, jj, 1, nt_a ) = zesrf 
     773 
     774            !++ Top Neumann B.C. 
     775            !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 ) 
     776            !z_elem_c ( ji,     jpka       ) = 0._wp 
     777            !z_elem_b ( ji,     jpka       ) = e3w_abl(jpka) - z_elem_a(ji, jpka ) 
     778            !tke_abl  ( ji, jj, jpka, nt_a ) = e3w_abl(jpka) * tke_abl( ji,jj, jpka, nt_n ) 
     779 
     780            !++ Top Dirichlet B.C. 
     781            z_elem_a ( ji,     jpka       ) = 0._wp 
     782            z_elem_c ( ji,     jpka       ) = 0._wp 
     783            z_elem_b ( ji,     jpka       ) = 1._wp 
     784            tke_abl  ( ji, jj, jpka, nt_a ) = zetop 
     785 
     786            zbn2 ( ji, jj,    1 ) = zbn2 ( ji, jj,      2 ) 
     787            zsh2 ( ji,        1 ) = zsh2 ( ji,          2 ) 
     788            zbn2 ( ji, jj, jpka ) = zbn2 ( ji, jj, jpkam1 ) 
     789            zsh2 ( ji,     jpka ) = zsh2 ( ji    , jpkam1 ) 
    771790         END DO 
    772791         !! 
     
    844863      ! Optional : could add pblh smoothing if pblh is noisy horizontally ... 
    845864      IF(ln_smth_pblh) THEN 
    846          CALL lbc_lnk( 'ablmod', pblh, 'T', 1., kfillmode = jpfillnothing) 
     865         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) !, kfillmode = jpfillnothing) 
    847866         CALL smooth_pblh( pblh, msk_abl ) 
    848          CALL lbc_lnk( 'ablmod', pblh, 'T', 1., kfillmode = jpfillnothing) 
     867         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.) !, kfillmode = jpfillnothing) 
    849868      ENDIF 
    850869      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    860879            ! 
    861880            DO ji = 1, jpi 
    862                mxld_abl ( ji, jj,    1 ) = mxl_min 
    863                mxld_abl ( ji, jj, jpka ) = mxl_min 
    864                mxlm_abl ( ji, jj,    1 ) = mxl_min 
    865                mxlm_abl ( ji, jj, jpka ) = mxl_min 
    866                zldw     ( ji,        1 ) = zrough(ji,jj) * rn_Lsfc 
    867                zlup     ( ji,     jpka ) = mxl_min 
     881               mxld_abl( ji, jj,    1 ) = mxl_min 
     882               mxld_abl( ji, jj, jpka ) = mxl_min 
     883               mxlm_abl( ji, jj,    1 ) = mxl_min 
     884               mxlm_abl( ji, jj, jpka ) = mxl_min 
     885               zldw    ( ji,        1 ) = zrough(ji,jj) * rn_Lsfc 
     886               zlup    ( ji,     jpka ) = mxl_min 
    868887            END DO 
    869888            ! 
     
    898917            DO jk = 1, jpka 
    899918               DO ji = 1, jpi 
     919!                  zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp) 
    900920                  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) 
    902921                  mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min ) 
    903                   !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ),  zlup( ji, jk ) ) 
    904922                  mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min ) 
    905923               END DO 
     
    971989#   undef zlup 
    972990#   undef zldw 
    973 !         ! 
     991         ! 
    974992      CASE ( 2 )           ! Bougeault & Lacarrere 89 length-scale 
    975993         ! 
     
    10071025            DO jkup=jk+1,jpka-1 
    10081026               DO ji = 1, jpi 
    1009                   zbuoy1             = MAX( zbn2(ji,jj,jkup  ), rsmall ) 
    1010                   zbuoy2             = MAX( zbn2(ji,jj,jkup-1), rsmall ) 
     1027                  zbuoy1 = MAX( zbn2(ji,jj,jkup  ), rsmall ) 
     1028                  zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall ) 
    10111029                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * & 
    10121030                     &               ( zbuoy1*(ghw_abl(jkup  )-ghw_abl(jk)) & 
     
    10381056            DO jkdwn=jk-1,1,-1 
    10391057               DO ji = 1, jpi 
    1040                   zbuoy1             = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 
    1041                   zbuoy2             = MAX( zbn2(ji,jj,jkdwn  ), rsmall ) 
     1058                  zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall ) 
     1059                  zbuoy2 = MAX( zbn2(ji,jj,jkdwn  ), rsmall ) 
    10421060                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  & 
    10431061                     &               * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) & 
     
    11661184                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) ) 
    11671185                     zldw(ji,jk) = zcff 
    1168                                          zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  ) 
     1186                     zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  ) 
    11691187                     ln_foundl(ji) = .true. 
    11701188                  END IF 
     
    11861204#   undef zlup 
    11871205#   undef zldw 
    1188          ! 
    11891206         ! 
    11901207         ! 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/par_abl.F90

    r12588 r12749  
    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     
     30   !LOGICAL , PUBLIC            ::   ln_topbc_neumann = .FALSE.  !: idealised testcases only 
    3131 
    3232   LOGICAL           , PUBLIC  ::   ln_rstart_abl    !: (de)activate abl restart 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ABL/sbcabl.F90

    r12592 r12749  
    168168      rn_Esfc =  1._wp / SQRT(rn_cm*rn_ceps) 
    169169      rn_Lsfc = vkarmn * SQRT(SQRT(rn_cm*rn_ceps)) / rn_cm 
    170       
     170 
    171171      IF(lwp) THEN 
    172172         WRITE(numout,*) 
     
    175175         IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale ' 
    176176         IF(nn_amxl==1) WRITE(numout,*) 'Modified Deardorff 80 length-scale ' 
    177          IF(nn_amxl==1) WRITE(numout,*) 'Bougeault and Lacarrere length-scale '       
     177         IF(nn_amxl==2) WRITE(numout,*) 'Bougeault and Lacarrere length-scale '       
    178178         IF(nn_amxl==3) WRITE(numout,*) 'Rodier et al. length-scale '    
    179179         WRITE(numout,*) ' Minimum value of atmospheric TKE           = ',tke_min,' m^2 s^-2' 
Note: See TracChangeset for help on using the changeset viewer.