Ignore:
Timestamp:
2019-12-12T09:30:08+01:00 (17 months ago)
Author:
laurent
Message:

Catches up with SBCBLK bugfixes of branch "dev_r12072_MERGE_OPTION2_2019"

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/ABL/sbcabl.F90

    r12193 r12199  
    22   !!====================================================================== 
    33   !!                       ***  MODULE  sbcabl  *** 
    4    !! Ocean forcing:  momentum, heat and freshwater flux formulation  
    5    !!                 derived from an ABL model  
     4   !! Ocean forcing:  momentum, heat and freshwater flux formulation 
     5   !!                 derived from an ABL model 
    66   !!===================================================================== 
    77   !! History :  4.0  !  2019-03  (F. Lemarié & G. Samson)  Original code 
     
    99 
    1010   !!---------------------------------------------------------------------- 
    11    !!   sbc_abl_init  : Initialization of ABL model based on namelist options  
    12    !!   sbc_abl       : driver for the computation of momentum, heat and freshwater  
     11   !!   sbc_abl_init  : Initialization of ABL model based on namelist options 
     12   !!   sbc_abl       : driver for the computation of momentum, heat and freshwater 
    1313   !!                   fluxes over ocean via the ABL model 
    1414   !!---------------------------------------------------------------------- 
     
    3535   USE ice    , ONLY : u_ice, v_ice, tm_su, ato_i      ! ato_i = total open water fractional area 
    3636   USE sbc_ice, ONLY : wndm_ice, utau_ice, vtau_ice 
    37 #endif    
     37#endif 
    3838#if ! defined key_iomput 
    3939   USE diawri    , ONLY : dia_wri_alloc_abl 
     
    5959      !! 
    6060      !! ** Purposes :   - read namelist section namsbc_abl 
    61       !!                 - initialize and check parameter values  
    62       !!                 - initialize variables of ABL model  
     61      !!                 - initialize and check parameter values 
     62      !!                 - initialize variables of ABL model 
    6363      !! 
    6464      !!---------------------------------------------------------------------- 
     
    7474         &                 rn_ldyn_min , rn_ldyn_max, rn_ltra_min, rn_ltra_max,   & 
    7575         &                 nn_amxl, rn_cm, rn_ct, rn_ce, rn_ceps, rn_Rod, rn_Ric, & 
    76          &                 ln_smth_pblh        
    77       !!--------------------------------------------------------------------- 
    78  
     76         &                 ln_smth_pblh 
     77      !!--------------------------------------------------------------------- 
     78 
     79      REWIND( numnam_ref )              ! Namelist namsbc_abl in reference namelist : ABL parameters 
    7980      READ  ( numnam_ref, namsbc_abl, IOSTAT = ios, ERR = 901 ) 
    8081901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in reference namelist' ) 
    8182      ! 
     83      REWIND( numnam_cfg )              ! Namelist namsbc_abl in configuration namelist : ABL parameters 
    8284      READ  ( numnam_cfg, namsbc_abl, IOSTAT = ios, ERR = 902 ) 
    8385902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in configuration namelist' ) 
     
    8587      IF(lwm) WRITE( numond, namsbc_abl ) 
    8688      ! 
    87       ! Check ABL mixing length option       
     89      ! Check ABL mixing length option 
    8890      IF( nn_amxl  < 0   .OR.  nn_amxl  > 2 )   & 
    89          &                 CALL ctl_stop( 'abl_init : bad flag, nn_amxl must be  0, 1 or 2 ' )          
    90       ! 
    91       ! Check ABL dyn restore option           
     91         &                 CALL ctl_stop( 'abl_init : bad flag, nn_amxl must be  0, 1 or 2 ' ) 
     92      ! 
     93      ! Check ABL dyn restore option 
    9294      IF( nn_dyn_restore  < 0   .OR.  nn_dyn_restore  > 2 )   & 
    93          &                 CALL ctl_stop( 'abl_init : bad flag, nn_dyn_restore must be  0, 1 or 2 ' )             
     95         &                 CALL ctl_stop( 'abl_init : bad flag, nn_dyn_restore must be  0, 1 or 2 ' ) 
    9496 
    9597      !!--------------------------------------------------------------------- 
    9698      !! Control prints 
    97       !!---------------------------------------------------------------------  
     99      !!--------------------------------------------------------------------- 
    98100      IF(lwp) THEN                              ! Control print (other namelist variable) 
    99101         WRITE(numout,*) 
     
    104106            IF(ln_geos_winds) THEN 
    105107               ln_geos_winds = .FALSE. 
    106                WRITE(numout,*) '      ABL -- geostrophic guide disabled (not compatible with ln_hpgls_frc = .T.)'                   
     108               WRITE(numout,*) '      ABL -- geostrophic guide disabled (not compatible with ln_hpgls_frc = .T.)' 
    107109            END IF 
    108110         ELSE IF( ln_geos_winds ) THEN 
    109             WRITE(numout,*) '      ABL -- winds forced by geostrophic winds'                
     111            WRITE(numout,*) '      ABL -- winds forced by geostrophic winds' 
    110112         ELSE 
    111             WRITE(numout,*) '      ABL -- Geostrophic winds and large-scale pressure gradient are ignored'                
     113            WRITE(numout,*) '      ABL -- Geostrophic winds and large-scale pressure gradient are ignored' 
    112114         END IF 
    113115         ! 
    114116         SELECT CASE ( nn_dyn_restore ) 
    115          CASE ( 0 )    
    116             WRITE(numout,*) '      ABL -- No restoring for ABL winds'    
     117         CASE ( 0 ) 
     118            WRITE(numout,*) '      ABL -- No restoring for ABL winds' 
    117119         CASE ( 1 ) 
    118             WRITE(numout,*) '      ABL -- Restoring of ABL winds only in the equatorial region '              
     120            WRITE(numout,*) '      ABL -- Restoring of ABL winds only in the equatorial region ' 
    119121         CASE ( 2 ) 
    120             WRITE(numout,*) '      ABL -- Restoring of ABL winds activated everywhere '               
     122            WRITE(numout,*) '      ABL -- Restoring of ABL winds activated everywhere ' 
    121123         END SELECT 
    122124         ! 
    123        IF( ln_smth_pblh )  WRITE(numout,*) '      ABL -- Smoothing of PBL height is activated'    
    124       ! 
     125         IF( ln_smth_pblh )  WRITE(numout,*) '      ABL -- Smoothing of PBL height is activated' 
     126        ! 
    125127      ENDIF 
    126128 
    127129      !!--------------------------------------------------------------------- 
    128130      !! Convert nudging coefficient from hours to 1/sec 
    129       !!---------------------------------------------------------------------          
     131      !!--------------------------------------------------------------------- 
    130132      zcff        = 1._wp / 3600._wp 
    131133      rn_ldyn_min = zcff / rn_ldyn_min 
    132       rn_ldyn_max = zcff / rn_ldyn_max       
     134      rn_ldyn_max = zcff / rn_ldyn_max 
    133135      rn_ltra_min = zcff / rn_ltra_min 
    134       rn_ltra_max = zcff / rn_ltra_max    
    135  
    136       !!--------------------------------------------------------------------- 
    137       !! ABL grid initialization  
    138       !!---------------------------------------------------------------------       
    139       CALL iom_open( TRIM(cn_dir)//TRIM(cn_dom), inum )  
     136      rn_ltra_max = zcff / rn_ltra_max 
     137 
     138      !!--------------------------------------------------------------------- 
     139      !! ABL grid initialization 
     140      !!--------------------------------------------------------------------- 
     141      CALL iom_open( TRIM(cn_dir)//TRIM(cn_dom), inum ) 
    140142      id     = iom_varid( inum, 'e3t_abl', kdimsz=idimsz, kndims=indims, lduld=lluldl ) 
    141       jpka   = idimsz(indims - COUNT( (/lluldl/) ) )  
     143      jpka   = idimsz(indims - COUNT( (/lluldl/) ) ) 
    142144      jpkam1 = jpka - 1 
    143145 
     
    150152 
    151153#if ! defined key_iomput 
    152      IF( dia_wri_alloc_abl()  /= 0 ) CALL ctl_stop( 'STOP', 'abl_init : unable to allocate arrays' ) 
     154      IF( dia_wri_alloc_abl()  /= 0 ) CALL ctl_stop( 'STOP', 'abl_init : unable to allocate arrays' ) 
    153155#endif 
    154156 
     
    156158         WRITE(numout,*) 
    157159         WRITE(numout,*) '    sbc_abl_init   : ABL Reference vertical grid' 
    158          WRITE(numout,*) '    ~~~~~~~'       
     160         WRITE(numout,*) '    ~~~~~~~' 
    159161         WRITE(numout, "(9x,'  level ght_abl   ghw_abl   e3t_abl   e3w_abl  ')" ) 
    160162         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, ght_abl(jk), ghw_abl(jk), e3t_abl(jk), e3w_abl(jk), jk = 1, jpka ) 
     
    162164 
    163165      !!--------------------------------------------------------------------- 
    164       !! Check TKE closure parameters  
    165       !!---------------------------------------------------------------------  
     166      !! Check TKE closure parameters 
     167      !!--------------------------------------------------------------------- 
    166168      rn_Sch  = rn_ce / rn_cm 
    167169      mxl_min = (avm_bak / rn_cm) / sqrt( tke_min ) 
     
    170172         WRITE(numout,*) 
    171173         WRITE(numout,*) '    abl_zdf_tke   : ABL TKE turbulent closure' 
    172          WRITE(numout,*) '    ~~~~~~~~~~~'       
     174         WRITE(numout,*) '    ~~~~~~~~~~~' 
    173175         IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale ' 
    174176         IF(nn_amxl==1) WRITE(numout,*) 'length-scale based on the distance to the PBL height ' 
    175177         WRITE(numout,*) ' Minimum value of atmospheric TKE           = ',tke_min,' m^2 s^-2' 
    176178         WRITE(numout,*) ' Minimum value of atmospheric mixing length = ',mxl_min,' m' 
    177          WRITE(numout,*) ' Constant for turbulent viscosity           = ',rn_Cm          
    178          WRITE(numout,*) ' Constant for turbulent diffusivity         = ',rn_Ct          
    179          WRITE(numout,*) ' Constant for Schmidt number                = ',rn_Sch             
    180          WRITE(numout,*) ' Constant for TKE dissipation               = ',rn_Ceps                         
     179         WRITE(numout,*) ' Constant for turbulent viscosity           = ',rn_Cm 
     180         WRITE(numout,*) ' Constant for turbulent diffusivity         = ',rn_Ct 
     181         WRITE(numout,*) ' Constant for Schmidt number                = ',rn_Sch 
     182         WRITE(numout,*) ' Constant for TKE dissipation               = ',rn_Ceps 
    181183      END IF 
    182184 
    183185      !!------------------------------------------------------------------------------------------- 
    184186      !! Compute parameters to build the vertical profile for the nudging term (used in abl_stp()) 
    185       !!-------------------------------------------------------------------------------------------        
     187      !!------------------------------------------------------------------------------------------- 
    186188      zcff1       = 1._wp / ( jp_bmax - jp_bmin )**3 
    187189      ! for active tracers 
     
    190192      jp_alp1_tra = -6._wp * zcff1 *  jp_bmax * jp_bmin  * ( rn_ltra_max - rn_ltra_min ) 
    191193      jp_alp0_tra = zcff1 * (  rn_ltra_max * jp_bmin*jp_bmin * (3._wp*jp_bmax - jp_bmin)      & 
    192          &                   - rn_ltra_min * jp_bmax*jp_bmax * (3._wp*jp_bmin - jp_bmax) )   
     194         &                   - rn_ltra_min * jp_bmax*jp_bmax * (3._wp*jp_bmin - jp_bmax) ) 
    193195      ! for dynamics 
    194196      jp_alp3_dyn = -2._wp * zcff1                       * ( rn_ldyn_max - rn_ldyn_min ) 
     
    196198      jp_alp1_dyn = -6._wp * zcff1 *  jp_bmax * jp_bmin  * ( rn_ldyn_max - rn_ldyn_min ) 
    197199      jp_alp0_dyn = zcff1 * (  rn_ldyn_max * jp_bmin*jp_bmin * (3._wp*jp_bmax - jp_bmin)      & 
    198          &                   - rn_ldyn_min * jp_bmax*jp_bmax * (3._wp*jp_bmin - jp_bmax) )        
     200         &                   - rn_ldyn_min * jp_bmax*jp_bmax * (3._wp*jp_bmin - jp_bmax) ) 
    199201 
    200202      jp_pblh_min = ghw_abl(     4) / jp_bmin  !<-- at least 3 grid points at the bottom have value rn_ltra_min 
     
    227229         &    + jp_alp1_tra * jp_bmin    + jp_alp0_tra ) * rdt_abl 
    228230      zcff1 = ( jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2   & 
    229          &    + jp_alp1_tra * jp_bmax    + jp_alp0_tra ) * rdt_abl  
     231         &    + jp_alp1_tra * jp_bmax    + jp_alp0_tra ) * rdt_abl 
    230232      IF(lwp) THEN 
    231233         WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff 
     
    242244      !!------------------------------------------------------------------------------------------- 
    243245      !! Initialize Coriolis frequency, equatorial restoring and land/sea mask 
    244       !!-------------------------------------------------------------------------------------------   
     246      !!------------------------------------------------------------------------------------------- 
    245247      fft_abl(:,:) = 2._wp * omega * SIN( rad * gphit(:,:) ) 
    246248 
    247249      ! Equatorial restoring 
    248250      IF( nn_dyn_restore == 1 ) THEN 
    249          zcff         = 2._wp * omega * SIN( rad * 90._wp )   !++ fmax  
     251         zcff         = 2._wp * omega * SIN( rad * 90._wp )   !++ fmax 
    250252         rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff) / zcff ) )**8 
    251253         !!GS: alternative shape 
     
    258260      msk_abl(:,:) = tmask(:,:,1) 
    259261 
    260       !!-------------------------------------------------------------------------------------------  
     262      !!------------------------------------------------------------------------------------------- 
    261263 
    262264      ! initialize 2D bulk fields AND 3D abl data 
     
    273275         CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step 
    274276 
    275           u_abl(:,:,:,nt_n      ) = sf(jp_wndi)%fnow(:,:,:) 
    276           v_abl(:,:,:,nt_n      ) = sf(jp_wndj)%fnow(:,:,:) 
     277         u_abl(:,:,:,nt_n      ) = sf(jp_wndi)%fnow(:,:,:) 
     278         v_abl(:,:,:,nt_n      ) = sf(jp_wndj)%fnow(:,:,:) 
    277279         tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:) 
    278280         tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:) 
    279     
     281 
    280282         tke_abl(:,:,:,nt_n     ) = tke_min 
    281283         avm_abl(:,:,:          ) = avm_bak 
    282284         avt_abl(:,:,:          ) = avt_bak 
    283          mxl_abl(:,:,:          ) = mxl_min  
    284          pblh   (:,:            ) = ghw_abl( 3 )  !<-- assume that the pbl contains 3 grid points  
     285         mxl_abl(:,:,:          ) = mxl_min 
     286         pblh   (:,:            ) = ghw_abl( 3 )  !<-- assume that the pbl contains 3 grid points 
    285287         u_abl  (:,:,:,nt_a     ) = 0._wp 
    286288         v_abl  (:,:,:,nt_a     ) = 0._wp 
     
    290292 
    291293      rhoa(:,:) = rho_air( tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), sf(jp_slp)%fnow(:,:,1) ) !!GS: rhoa must be (re)computed here here to avoid division by zero in blk_ice_1 (TBI) 
    292       
     294 
    293295   END SUBROUTINE sbc_abl_init 
    294296 
     
    299301      !! 
    300302      !! ** Purpose :   provide the momentum, heat and freshwater fluxes at 
    301       !!      the ocean surface from an ABL calculation at each oceanic time step   
    302       !! 
    303       !! ** Method  :  
     303      !!      the ocean surface from an ABL calculation at each oceanic time step 
     304      !! 
     305      !! ** Method  : 
    304306      !!              - Pre-compute part of turbulent fluxes in blk_oce_1 
    305       !!              - Perform 1 time-step of the ABL model    
     307      !!              - Perform 1 time-step of the ABL model 
    306308      !!              - Finalize flux computation in blk_oce_2 
    307309      !! 
     
    320322#if defined key_si3 
    321323      REAL(wp), DIMENSION(jpi,jpj) ::   zssqi, zcd_dui, zseni, zevpi 
    322 #endif      
     324#endif 
    323325      INTEGER                      ::   jbak, jbak_dta, ji, jj 
    324326      !!--------------------------------------------------------------------- 
     
    327329      !! 1 - Read Atmospheric 3D data for large-scale forcing 
    328330      !!------------------------------------------------------------------------------------------- 
    329        
     331 
    330332      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step 
    331       
     333 
    332334      !!------------------------------------------------------------------------------------------- 
    333335      !! 2 - Compute Cd x ||U||, Ch x ||U||, Ce x ||U||, and SSQ using now fields 
    334       !!-------------------------------------------------------------------------------------------   
     336      !!------------------------------------------------------------------------------------------- 
    335337 
    336338      CALL blk_oce_1( kt,  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),   &   !   <<= in 
     
    338340         &                sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m     ,   &   !   <<= in 
    339341         &                sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) ,   &   !   <<= in 
    340          &                zssq, zcd_du, zsen, zevp                          )       !   =>> out 
    341    
     342         &                tsk_m, zssq, zcd_du, zsen, zevp                       )   !   =>> out 
     343 
    342344#if defined key_si3 
    343345      CALL blk_ice_1(  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),    &   !   <<= in 
     
    345347         &            sf(jp_slp)%fnow(:,:,1)  ,  u_ice, v_ice, tm_su    ,    &   !   <<= in 
    346348         &            pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui )   !   <<= out 
    347 #endif   
     349#endif 
    348350 
    349351      !!------------------------------------------------------------------------------------------- 
    350352      !! 3 - Advance ABL variables from now (n) to after (n+1) 
    351       !!-------------------------------------------------------------------------------------------     
    352    
    353       CALL abl_stp( kt, sst_m, ssu_m, ssv_m, zssq,                          &   !   <<= in 
     353      !!------------------------------------------------------------------------------------------- 
     354 
     355      CALL abl_stp( kt, tsk_m, ssu_m, ssv_m, zssq,                          &   !   <<= in 
    354356         &              sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:),   &   !   <<= in 
    355357         &              sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:),   &   !   <<= in 
     
    358360         &              zcd_du, zsen, zevp,                                 &   !   <=> in/out 
    359361         &              wndm, utau, vtau, taum                              &   !   =>> out 
    360 #if defined key_si3           
     362#if defined key_si3 
    361363         &            , tm_su, u_ice, v_ice, zssqi, zcd_dui                 &   !   <<= in 
    362364         &            , zseni, zevpi, wndm_ice, ato_i                       &   !   <<= in 
    363365         &            , utau_ice, vtau_ice                                  &   !   =>> out 
    364 #endif            
     366#endif 
    365367         &                                                                  ) 
    366368      !!------------------------------------------------------------------------------------------- 
    367       !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since  
     369      !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since 
    368370      !!                                                                time swap is done in abl_stp 
    369       !!-------------------------------------------------------------------------------------------    
     371      !!------------------------------------------------------------------------------------------- 
    370372 
    371373      CALL blk_oce_2( tq_abl(:,:,2,nt_n,jp_ta),                            & 
    372374         &            sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1),   & 
    373375         &            sf(jp_prec)%fnow(:,:,1) , sf(jp_snow)%fnow(:,:,1),   & 
    374          &            sst_m, zsen, zevp                                ) 
    375  
    376       CALL abl_rst_opn( kt )                       ! Open abl restart file (if necessary)  
    377       IF( lrst_abl ) CALL abl_rst_write( kt )      ! -- abl restart file  
     376         &            tsk_m, zsen, zevp                                ) 
     377 
     378      CALL abl_rst_opn( kt )                       ! Open abl restart file (if necessary) 
     379      IF( lrst_abl ) CALL abl_rst_write( kt )      ! -- abl restart file 
    378380 
    379381#if defined key_si3 
    380382      ! Avoid a USE abl in icesbc module 
    381383      sf(jp_tair)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_ta);  sf(jp_humi)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_qa) 
    382 #endif  
     384#endif 
    383385 
    384386   END SUBROUTINE sbc_abl 
Note: See TracChangeset for help on using the changeset viewer.