Changeset 11363


Ignore:
Timestamp:
2019-07-29T16:47:01+02:00 (13 months ago)
Author:
gsamson
Message:

dev_r11265_ABL : see #2131

  • fix xml files to get xios working with ABL
  • introduce rdt_abl = rdt * nn_fsbc to get ABL working with nn_fsbc > 1
  • does not change results compared to rev11360 with nn_fsbc = 1
Location:
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-ice.xml

    r10911 r11363  
    9090          <field id="qemp_ice"     long_name="Downward Heat Flux from E-P over ice"                                                                              unit="W/m2" /> 
    9191          <field id="albedo"       long_name="Mean albedo over sea ice and ocean"                                                                                unit=""     /> 
     92          <field id="Cd_ice"       long_name="Momentum turbulent exchange coefficient"                                                                           unit=""     /> 
     93          <field id="Ch_ice"       long_name="Heat turbulent exchange coefficient"                                                                               unit=""     /> 
    9294      
    9395     <!-- trends --> 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-oce.xml

    r11306 r11363  
    2424 
    2525        <field id="toce"         long_name="temperature"                         standard_name="sea_water_potential_temperature"   unit="degC"     grid_ref="grid_T_3D"/> 
    26         <field id="toce_e3t"     long_name="temperature (thickness weighted)"                                                      unit="degC"     grid_ref="grid_T_3D" > toce * e3t </field > 
     26        <field id="toce_e3t"     long_name="temperature (thickness weighted)"                                                      unit="degC"     grid_ref="grid_T_3D" > toce * e3t </field> 
    2727        <field id="soce"         long_name="salinity"                            standard_name="sea_water_practical_salinity"      unit="1e-3"     grid_ref="grid_T_3D"/> 
    28         <field id="soce_e3t"     long_name="salinity    (thickness weighted)"                                                      unit="1e-3"     grid_ref="grid_T_3D" > soce * e3t </field > 
     28        <field id="soce_e3t"     long_name="salinity    (thickness weighted)"                                                      unit="1e-3"     grid_ref="grid_T_3D" > soce * e3t </field> 
    2929 
    3030        <!-- t-eddy viscosity coefficients (ldfdyn) --> 
     
    3333 
    3434        <field id="sst"          long_name="sea surface temperature"                            standard_name="sea_surface_temperature"             unit="degC"     /> 
    35         <field id="sst2"         long_name="square of sea surface temperature"                  standard_name="square_of_sea_surface_temperature"   unit="degC2"     > sst * sst </field > 
     35        <field id="sst2"         long_name="square of sea surface temperature"                  standard_name="square_of_sea_surface_temperature"   unit="degC2"     > sst * sst </field> 
    3636        <field id="sstmax"       long_name="max of sea surface temperature"   field_ref="sst"   operation="maximum"                                                 /> 
    3737        <field id="sstmin"       long_name="min of sea surface temperature"   field_ref="sst"   operation="minimum"                                                 /> 
     
    4545         
    4646        <field id="sss"          long_name="sea surface salinity"                               standard_name="sea_surface_salinity"                unit="1e-3"     /> 
    47         <field id="sss2"         long_name="square of sea surface salinity"                                                                         unit="1e-6"      > sss * sss </field > 
     47        <field id="sss2"         long_name="square of sea surface salinity"                                                                         unit="1e-6"      > sss * sss </field> 
    4848        <field id="sssmax"       long_name="max of sea surface salinity"      field_ref="sss"   operation="maximum"                                                 /> 
    4949        <field id="sssmin"       long_name="min of sea surface salinity"      field_ref="sss"   operation="minimum"                                                 /> 
     
    5454 
    5555        <field id="ssh"          long_name="sea surface height"                                 standard_name="sea_surface_height_above_geoid"             unit="m" /> 
    56         <field id="ssh2"         long_name="square of sea surface height"                       standard_name="square_of_sea_surface_height_above_geoid"   unit="m2" > ssh * ssh </field > 
     56        <field id="ssh2"         long_name="square of sea surface height"                       standard_name="square_of_sea_surface_height_above_geoid"   unit="m2" > ssh * ssh </field> 
    5757        <field id="wetdep"       long_name="wet depth"                                          standard_name="wet_depth"                                  unit="m" /> 
    5858        <field id="sshmax"       long_name="max of sea surface height"        field_ref="ssh"   operation="maximum"                                                 /> 
     
    352352      <field_group id="ABL" > <!-- time step automaticaly defined based on nn_fsbc --> 
    353353   <!-- variables available with ABL on atmospheric T grid--> 
    354    <field_group id="grid_ABL" grid_ref="grid_TA_3D" > 
     354   <field_group id="grid_ABL3D" grid_ref="grid_TA_3D" > 
    355355          <field id="u_abl"      long_name="ABL i-horizontal velocity"     standard_name="abl_x_velocity" unit="m/s"      /> 
    356356          <field id="v_abl"      long_name="ABL j-horizontal velocity"     standard_name="abl_y_velocity" unit="m/s"      /> 
    357           <field id="t_abl"      long_name="ABL potential temperature"  standard_name="abl_theta" unit="K"        /> 
    358           <field id="q_abl"      long_name="ABL specific humidity"         standard_name="abl_qspe" unit="kg/kg"    /> 
    359           <field id="pblh"       long_name="ABL height"                          standard_name="abl_height" unit="m"    grid_ref="grid_T_2D"     /> 
     357          <field id="t_abl"      long_name="ABL potential temperature"     standard_name="abl_theta"      unit="K"        /> 
     358          <field id="q_abl"      long_name="ABL specific humidity"         standard_name="abl_qspe"       unit="kg/kg"    /> 
     359          <field id="coeft"      long_name="ABL nudging coefficient"       standard_name="coeft"          unit=""         /> 
     360   </field_group> 
     361   <field_group id="grid_ABL2D" grid_ref="grid_TA_2D" > 
     362          <field id="pblh"       long_name="ABL height"                    standard_name="abl_height"     unit="m"        /> 
    360363   </field_group> 
    361364      </field_group> <!-- ABL --> 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/grid_def_nemo.xml

    r11275 r11363  
    5959 
    6060       <!-- ABL grid definition --> 
     61       <grid id="grid_TA_2D"> 
     62         <domain domain_ref="grid_T" /> 
     63       </grid> 
    6164       <grid id="grid_TA_3D"> 
    6265         <domain domain_ref="grid_T" /> 
    6366         <axis  id="ght_abl" /> 
    64        </grid> 
    65        <grid id="grid_TA_2D"> 
    66          <domain domain_ref="grid_T" /> 
    6767       </grid> 
    6868       <grid id="grid_WA_3D"> 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ref

    r11334 r11363  
    283283&namsbc_abl    !   Atmospheric Boundary Layer formulation           (ln_abl = T) 
    284284!----------------------------------------------------------------------- 
    285    cn_dir         = 'atm_forcing/'      !  root directory for the location of the ABL grid file 
     285   cn_dir         = './'      !  root directory for the location of the ABL grid file 
    286286   cn_dom         = 'dom_cfg_abl.nc' 
    287287   ln_hpgls_frc   = .false. 
     
    301301   rn_Rod        = 0.15       ! c0 in RMCA17 mixing length formulation (not yet implemented) 
    302302   rn_Ric        = 0.139      !  Critical Richardson number (to compute PBL height and diffusivities) 
    303    ln_smth_pblh  = .true.    !  Smoothing of PBL height with a 2d Hanning filter 
     303   ln_smth_pblh  = .false.   !  Smoothing of PBL height with a 2d Hanning filter 
    304304/ 
    305305!----------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90

    r11348 r11363  
    153153         DO jk = 3, jpkam1 
    154154            DO ji = 1, jpi   ! vector opt. 
    155                z_elem_a( ji,     jk              ) = - rdt * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )          ! lower-diagonal 
    156                z_elem_c( ji,     jk              ) = - rdt * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )          ! upper-diagonal        
     155               z_elem_a( ji,     jk              ) = - rdt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )          ! lower-diagonal 
     156               z_elem_c( ji,     jk              ) = - rdt_abl * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )          ! upper-diagonal        
    157157               z_elem_b( ji,     jk              ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )  !       diagonal            
    158158            END DO 
     
    162162            ! Neumann at the bottom           
    163163            z_elem_a( ji,     2              ) = 0._wp 
    164             z_elem_c( ji,     2              ) = - rdt * Avt_abl( ji, jj, 2   ) / e3w_abl( 2   )                                          
     164            z_elem_c( ji,     2              ) = - rdt_abl * Avt_abl( ji, jj, 2   ) / e3w_abl( 2   )                                          
    165165            ! Homogeneous Neumann at the top 
    166             z_elem_a( ji,     jpka           ) = - rdt * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )  
     166            z_elem_a( ji,     jpka           ) = - rdt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )  
    167167            z_elem_c( ji,     jpka           ) = 0._wp 
    168168            z_elem_b( ji,     jpka           ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka ) 
     
    185185                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj) 
    186186#endif               
    187               z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1             
    188                   tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt * zztmp2                
     187              z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1               
     188                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt_abl * zztmp2                
    189189              tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra ) 
    190190               END DO  
     
    197197              zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)     
    198198#endif    
    199                   z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1 
    200                   tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt * zztmp2                
     199                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1 
     200                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt_abl * zztmp2                
    201201                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra ) 
    202202               END DO                                      
     
    244244         DO jj = 1, jpj 
    245245            DO ji = 1, jpi            
    246                zcff = ( fft_abl(ji,jj) * rdt )*( fft_abl(ji,jj) * rdt )  ! (f dt)**2 
     246               zcff = ( fft_abl(ji,jj) * rdt_abl )*( fft_abl(ji,jj) * rdt_abl )  ! (f dt)**2 
    247247       
    248248               u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(  & 
    249249                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n )    & 
    250                   &                 +  rdt * fft_abl(ji, jj) * v_abl ( ji , jj  , jk, nt_n ) )  & 
     250                  &                 +  rdt_abl * fft_abl(ji, jj) * v_abl ( ji , jj  , jk, nt_n ) )  & 
    251251                  &                               / (1._wp + gamma_Cor*gamma_Cor*zcff) 
    252252                   
    253253               v_abl( ji, jj, jk, nt_a ) =  e3t_abl(jk) *(  & 
    254254                  &        (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n )   & 
    255                   &                 -  rdt * fft_abl(ji, jj) * u_abl ( ji   , jj, jk, nt_n )  ) & 
     255                  &                 -  rdt_abl * fft_abl(ji, jj) * u_abl ( ji   , jj, jk, nt_n )  ) & 
    256256                  &                                / (1._wp + gamma_Cor*gamma_Cor*zcff)                 
    257257            END DO 
     
    267267               DO ji = 1, jpi  
    268268                  u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   & 
    269                      &                      - rdt * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk) 
     269                     &                      - rdt_abl * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk) 
    270270                  v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   & 
    271                      &                      + rdt * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk) 
     271                     &                      + rdt_abl * e3t_abl(jk) * fft_abl(ji, jj  ) * pgu_dta(ji  ,jj  ,jk) 
    272272               END DO 
    273273            END DO 
     
    280280            DO jk = 1, jpka 
    281281               DO ji = 1, jpi  
    282                   u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgu_dta(ji,jj,jk) 
    283                   v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgv_dta(ji,jj,jk)   
     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)   
    284284               ENDDO 
    285285            ENDDO 
     
    298298         DO jk = 3, jpkam1 
    299299            DO ji = 1, jpi   
    300                z_elem_a( ji,     jk ) = - rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal 
    301                z_elem_c( ji,     jk ) = - rdt * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal                 
     300               z_elem_a( ji,     jk ) = - rdt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )  ! lower-diagonal 
     301               z_elem_c( ji,     jk ) = - rdt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )  ! upper-diagonal                 
    302302               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                             !       diagonal 
    303303            END DO 
     
    307307            !++ Surface boundary condition 
    308308            z_elem_a( ji,     2    ) = 0._wp 
    309             z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                        
     309            z_elem_c( ji,     2    ) = - rdt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                        
    310310            ! 
    311311         zztmp1  = pcd_du(ji, jj) 
     
    316316         zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 
    317317#endif            
    318          z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1           
    319             u_abl( ji, jj,    2, nt_a ) =      u_abl( ji, jj,    2, nt_a ) + rdt * zztmp2 
     318         z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1          
     319            u_abl( ji, jj,    2, nt_a ) =      u_abl( ji, jj,    2, nt_a ) + rdt_abl * zztmp2 
    320320          
    321321            !++ Top Neumann B.C. 
    322             !z_elem_a( ji,     jpka ) = - 0.5_wp * rdt * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka )  
     322            !z_elem_a( ji,     jpka ) = - 0.5_wp * rdt_abl * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka )  
    323323            !z_elem_c( ji,     jpka ) = 0._wp 
    324324            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                                
     
    365365         DO jk = 3, jpkam1 
    366366            DO ji = 1, jpi    
    367                z_elem_a( ji,     jk ) = -rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal 
    368                z_elem_c( ji,     jk ) = -rdt * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal               
     367               z_elem_a( ji,     jk ) = -rdt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal 
     368               z_elem_c( ji,     jk ) = -rdt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal               
    369369               z_elem_b( ji,     jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )                              !       diagonal 
    370370            END DO 
     
    374374            !++ Surface boundary condition 
    375375            z_elem_a( ji,     2    ) = 0._wp 
    376             z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )         
     376            z_elem_c( ji,     2    ) = - rdt_abl * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )         
    377377            ! 
    378378         zztmp1 = pcd_du(ji, jj) 
     
    383383         zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 
    384384#endif          
    385             z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1  
    386             v_abl( ji, jj,    2, nt_a ) =         v_abl( ji, jj, 2, nt_a ) + rdt * zztmp2 
     385            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1  
     386            v_abl( ji, jj,    2, nt_a ) =         v_abl( ji, jj, 2, nt_a ) + rdt_abl * zztmp2 
    387387            !++ Top Neumann B.C.             
    388             !z_elem_a( ji,     jpka ) = -rdt * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )  
     388            !z_elem_a( ji,     jpka ) = -rdt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )  
    389389            !z_elem_c( ji,     jpka ) = 0._wp 
    390390            !z_elem_b( ji,     jpka ) = e3t_abl( jpka ) - z_elem_a( ji,     jpka )                                                
     
    434434            DO jj = 2, jpj 
    435435               DO ji = 2, jpi 
    436                   zcff1  = pblh( ji, jj ) 
    437                   zsig   = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )                         
    438                   zsig  =            MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )  
     436                  zcff1 = pblh( ji, jj ) 
     437                  zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )                         
     438                  zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )  
    439439                  zmsk  = msk_abl(ji,jj) 
    440440                  zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2   & 
    441441                     &  + jp_alp1_dyn * zsig    + jp_alp0_dyn 
    442                   zcff  = (1._wp-zmsk) + zmsk * rdt * zcff2   ! zcff = 1 for masked points 
    443                    
    444               zcff  = zcff * rest_eq(ji,jj) ; z_cft( ji, jj, jk ) = zcff 
     442                  zcff  = (1._wp-zmsk) + zmsk * zcff2 * rdt   ! zcff = 1 for masked points 
     443                                                              ! rdt = rdt_abl / nn_fsbc                           
     444        zcff  = zcff * rest_eq(ji,jj) ; z_cft( ji, jj, jk ) = zcff 
    445445               
    446446                  u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  u_abl( ji, jj, jk, nt_a )   & 
     
    448448                  v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  v_abl( ji, jj, jk, nt_a )   & 
    449449                     &                               + zcff   * pv_dta( ji, jj, jk       ) 
    450                 END DO 
     450               END DO 
    451451            END DO    
    452452         !------------- 
     
    460460         DO jj = 1,jpj 
    461461            DO ji = 1,jpi  
    462                zcff1  = pblh( ji, jj ) 
    463                zsig   = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) ) 
    464                zsig   =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )  
    465                zmsk   = msk_abl(ji,jj) 
    466                zcff2  = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2   & 
    467                   &   + jp_alp1_tra * zsig    + jp_alp0_tra 
    468                zcff   = (1._wp-zmsk) + zmsk * rdt * zcff2   ! zcff = 1 for masked points       
     462               zcff1 = pblh( ji, jj ) 
     463               zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) ) 
     464               zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )  
     465               zmsk  = msk_abl(ji,jj) 
     466               zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2   & 
     467                  &  + jp_alp1_tra * zsig    + jp_alp0_tra 
     468               zcff  = (1._wp-zmsk) + zmsk * zcff2 * rdt   ! zcff = 1 for masked points 
     469                                                           ! rdt = rdt_abl / nn_fsbc                           
    469470               tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta )   & 
    470471                  &                                       + zcff   * pt_dta( ji, jj, jk ) 
     
    483484      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    484485     ! 
    485       CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a), 'T', -1., v_abl(:,:,:,nt_a), 'T', -1., kfillmode = jpfillnothing  ) 
    486       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... 
     486      CALL lbc_lnk_multi( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1.,  v_abl(:,:,:,nt_a      ), 'T', -1., kfillmode = jpfillnothing  ) 
     487      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... 
    487488     ! 
    488489      CALL iom_put ( 'u_abl',  u_abl(:,:,2:jpka,nt_a      )) 
     
    490491      CALL iom_put ( 't_abl', tq_abl(:,:,2:jpka,nt_a,jp_ta)) 
    491492      CALL iom_put ( 'q_abl', tq_abl(:,:,2:jpka,nt_a,jp_qa)) 
    492 !      CALL iom_put ( 'coeft',  z_cft(:,:,2:jpka           )) 
     493      CALL iom_put ( 'coeft',  z_cft(:,:,2:jpka           )) 
    493494      CALL iom_put (  'pblh',   pblh(:,:                  )) 
    494495      ! 
     
    676677               zbuoy        = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )  
    677678                
    678                z_elem_a( ji,     jk )   = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk   ) ! lower-diagonal 
    679                z_elem_c( ji,     jk )   = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk   )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal           
     679               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 
     680               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           
    680681               IF( (zbuoy + zshear) .gt. 0.) THEN    ! Patankar trick to avoid negative values of TKE 
    681682                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
    682                      &                     + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)     ! diagonal        
    683                   tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt * ( zbuoy + zshear ) )             ! right-hand-side 
     683                     &                     + e3w_abl(jk) * rdt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)     ! diagonal        
     684                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt_abl * ( zbuoy + zshear ) )             ! right-hand-side 
    684685               ELSE 
    685686                  z_elem_b( ji,     jk )   = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   & 
    686                      &                     + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)   &  ! diagonal     
    687                      &                     - e3w_abl(jk) * rdt * zbuoy    
    688                   tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt *  zshear )             ! right-hand-side                      
     687                     &                     + e3w_abl(jk) * rdt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk)   &  ! diagonal     
     688                     &                     - e3w_abl(jk) * rdt_abl * zbuoy    
     689                  tke_abl( ji, jj, jk, nt_a )  = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt_abl *  zshear )             ! right-hand-side                      
    689690               END IF 
    690691            END DO 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/par_abl.F90

    r11322 r11363  
    1919   INTEGER , PUBLIC, PARAMETER ::   jp_qa  = 2     !: indice for humidity 
    2020   INTEGER , PUBLIC, PARAMETER ::   jptime = 2     !: number of time indices stored in memory 
     21 
    2122   !!--------------------------------------------------------------------- 
    2223   !! namABL Namelist options   
     
    6970   ! parameter for the semi-implicit treatment of Coriolis term   
    7071   REAL(wp), PUBLIC, PARAMETER ::   gamma_Cor  = 0.55_wp 
     72   ! ABL timestep 
     73   REAL(wp), PUBLIC            :: rdt_abl 
    7174 
    7275   !!--------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90

    r11360 r11363  
    200200      jp_pblh_max = ghw_abl(jpka-3) / jp_bmax  !<-- at least 3 grid points at the top    have value rn_ltra_max 
    201201 
     202      ! ABL timestep 
     203      rdt_abl = nn_fsbc * rdt 
     204 
    202205      ! Check parameters for dynamics 
    203       zcff  = jp_alp3_dyn * jp_bmin**3 + jp_alp2_dyn * jp_bmin**2   & 
    204          &  + jp_alp1_dyn * jp_bmin    + jp_alp0_dyn  
    205       zcff1 = jp_alp3_dyn * jp_bmax**3 + jp_alp2_dyn * jp_bmax**2   & 
    206          &  + jp_alp1_dyn * jp_bmax    + jp_alp0_dyn       
     206      zcff  = ( jp_alp3_dyn * jp_bmin**3 + jp_alp2_dyn * jp_bmin**2   & 
     207         &    + jp_alp1_dyn * jp_bmin    + jp_alp0_dyn ) * rdt_abl 
     208      zcff1 = ( jp_alp3_dyn * jp_bmax**3 + jp_alp2_dyn * jp_bmax**2   & 
     209         &    + jp_alp1_dyn * jp_bmax    + jp_alp0_dyn ) * rdt_abl 
    207210      IF(lwp) THEN 
    208211         IF(nn_dyn_restore > 0) THEN 
    209             WRITE(numout,*) ' ABL Minimum value for dynamics restoring = ',zcff  * rdt 
    210             WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1 * rdt 
     212            WRITE(numout,*) ' ABL Minimum value for dynamics restoring = ',zcff 
     213            WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1 
    211214            ! Check that restoring coefficients are between 0 and 1 
    212             IF( zcff1 * rdt > 1._wp .OR. zcff1 * rdt < 0._wp )   & 
     215            !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp )   & 
     216            IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp )   & 
    213217               &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_max' ) 
    214             IF( zcff  * rdt > 1._wp .OR. zcff  * rdt < 0._wp )   & 
     218            !IF( zcff  > 1._wp .OR. zcff  < 0._wp )   & 
     219            IF( zcff  > nn_fsbc .OR. zcff  < 0._wp )   & 
    215220               &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_min' ) 
    216             IF( zcff  * rdt > zcff1 * rdt ) & 
     221            IF( zcff  > zcff1                    )  & 
    217222               &                   CALL ctl_stop( 'abl_init : rn_ldyn_max should be smaller than rn_ldyn_min' ) 
    218223         END IF 
     
    220225 
    221226      ! Check parameters for active tracers 
    222       zcff  = jp_alp3_tra * jp_bmin**3 + jp_alp2_tra * jp_bmin**2   & 
    223          &  + jp_alp1_tra * jp_bmin    + jp_alp0_tra  
    224       zcff1 = jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2   & 
    225          &  + jp_alp1_tra * jp_bmax    + jp_alp0_tra  
     227      zcff  = ( jp_alp3_tra * jp_bmin**3 + jp_alp2_tra * jp_bmin**2   & 
     228         &    + jp_alp1_tra * jp_bmin    + jp_alp0_tra ) * rdt_abl 
     229      zcff1 = ( jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2   & 
     230         &    + jp_alp1_tra * jp_bmax    + jp_alp0_tra ) * rdt_abl  
    226231      IF(lwp) THEN 
    227          WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff  * rdt 
    228          WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1 * rdt 
     232         WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff 
     233         WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1 
    229234         ! Check that restoring coefficients are between 0 and 1 
    230          IF( zcff1 * rdt > 1._wp .OR. zcff1 * rdt < 0._wp )   & 
     235         !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp )   & 
     236         IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp )   & 
    231237            &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_max' ) 
    232          IF( zcff  * rdt > 1._wp .OR. zcff  * rdt < 0._wp )   & 
     238         !IF( zcff  > 1._wp .OR. zcff  < 0._wp )   & 
     239         IF( zcff  > nn_fsbc .OR. zcff  < 0._wp )   & 
    233240            &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_min' ) 
    234          IF( zcff  * rdt > zcff1  * rdt ) & 
     241         IF( zcff  > zcff1                    )  & 
    235242            &                   CALL ctl_stop( 'abl_init : rn_ltra_max should be smaller than rn_ltra_min' ) 
    236243      END IF 
     
    244251      IF( nn_dyn_restore == 1 ) THEN 
    245252         zcff         = 2._wp * omega * SIN( rad * 90._wp )   !++ fmax  
    246          rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff)/zcff ) )**8 
     253         rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff) / zcff ) )**8 
     254         !!GS: alternative shape 
     255         !rest_eq(:,:) = SIN( 0.5_wp*rpi*(zcff - ABS(ff_t(:,:))) / (zcff - 3.e-5) )**8 
     256         !WHERE(ABS(ff_t(:,:)).LE.3.e-5) rest_eq(:,:) = 1._wp 
    247257      ELSE 
    248258         rest_eq(:,:) = 1._wp 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90

    r11360 r11363  
    654654         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm 
    655655            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    656 ! 
     656         ! 
    657657         IF( ln_abl ) THEN 
    658          !                                                                                      !!! nid_A : 3D 
    659          CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl 
     658            CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl 
    660659               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
    661660            CALL histdef( nid_A, "q_abl", "Humidity"                  , "kg/kg"    ,       &  ! q_abl 
     
    677676               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout ) 
    678677#endif 
    679           CALL histend( nid_A, snc4chunks=snc4set ) 
    680        ! 
    681        ENDIF 
    682 ! 
     678            CALL histend( nid_A, snc4chunks=snc4set ) 
     679         ENDIF 
     680         ! 
    683681         IF( ln_icebergs ) THEN 
    684682            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , & 
     
    838836      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction    
    839837      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed    
    840 ! 
     838      ! 
    841839      IF( ln_abl ) THEN  
    842         ALLOCATE( zw3d_abl(jpi,jpj,jpka) ) 
    843         IF( ln_mskland )   THEN  
    844           DO jk=1,jpka 
    845              zw3d_abl(:,:,jk) = tmask(:,:,1) 
    846             END DO 
    847       ELSE 
     840         ALLOCATE( zw3d_abl(jpi,jpj,jpka) ) 
     841         IF( ln_mskland )   THEN  
     842            DO jk=1,jpka 
     843               zw3d_abl(:,:,jk) = tmask(:,:,1) 
     844            END DO        
     845        ELSE 
    848846            zw3d_abl(:,:,:) = 1._wp      
    849847         ENDIF        
    850       CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh  
    851         CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl 
    852         CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl 
    853         CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl 
    854         CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl      
    855         CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl 
    856         CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl 
    857         CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl   
     848        CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh  
     849         CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl 
     850         CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl 
     851         CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl 
     852         CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl        
     853         CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl 
     854         CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl 
     855         CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl  
    858856#if defined key_si3 
    859857         CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i 
    860858#endif 
    861       DEALLOCATE(zw3d_abl) 
    862      ENDIF 
    863 ! 
     859        DEALLOCATE(zw3d_abl) 
     860      ENDIF 
     861      ! 
    864862      IF( ln_icebergs ) THEN 
    865863         ! 
Note: See TracChangeset for help on using the changeset viewer.