Changeset 12199


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

Catches up with SBCBLK bugfixes of branch "dev_r12072_MERGE_OPTION2_2019"

Location:
NEMO/branches/2019/dev_r11943_MERGE_2019
Files:
1 deleted
11 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 
  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/abl.F90

    r12182 r12199  
    22   !!====================================================================== 
    33   !!                      ***  MODULE  abl  *** 
    4    !! Abl        :  ABL dynamics and active tracers defined in memory  
     4   !! Abl        :  ABL dynamics and active tracers defined in memory 
    55   !!====================================================================== 
    66   USE par_kind        ! abl parameters 
    7    
     7 
    88   IMPLICIT NONE 
    99   PRIVATE 
    10    !! --------------------------                            !  
     10   !! --------------------------                            ! 
    1111   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:)   ::   u_abl        !: i-horizontal velocity   [m/s] 
    1212   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:)   ::   v_abl        !: j-horizontal velocity   [m/s] 
    13    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:,:) ::   tq_abl       !: 4D T-q fields           [Kelvin,kg/kg]  
     13   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:,:) ::   tq_abl       !: 4D T-q fields           [Kelvin,kg/kg] 
    1414   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   avm_abl      !: turbulent viscosity   [m2/s] 
    1515   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     ::   avt_abl      !: turbulent diffusivity [m2/s] 
     
    2121   ! 
    2222   INTEGER , PUBLIC :: nt_n, nt_a       !: now / after indices (equal 1 or 2) 
    23    !  
     23   ! 
    2424   !!---------------------------------------------------------------------- 
    2525   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
    26    !! $Id: abl.F90 4990 2014-12-15 16:42:49Z timgraham $  
     26   !! $Id: abl.F90 4990 2014-12-15 16:42:49Z timgraham $ 
    2727   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    2828   !!---------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/sbc_oce.F90

    r12182 r12199  
    22   !!====================================================================== 
    33   !!                       ***  MODULE  sbc_oce  *** 
    4    !! Surface module :   variables defined in core memory  
     4   !! Surface module :   variables defined in core memory 
    55   !!====================================================================== 
    66   !! History :  3.0  ! 2006-06  (G. Madec)  Original code 
     
    99   !!             -   ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step 
    1010   !!            3.3  ! 2010-10  (J. Chanut, C. Bricaud)  add the surface pressure forcing 
    11    !!            4.0  ! 2012-05  (C. Rousset) add attenuation coef for use in ice model  
     11   !!            4.0  ! 2012-05  (C. Rousset) add attenuation coef for use in ice model 
    1212   !!            4.0  ! 2016-06  (L. Brodeau) new unified bulk routine (based on AeroBulk) 
    13    !!            4.0  ! 2019-03  (F. Lemarié, G. Samson) add compatibility with ABL mode     
     13   !!            4.0  ! 2019-03  (F. Lemarié, G. Samson) add compatibility with ABL mode 
    1414   !!---------------------------------------------------------------------- 
    1515 
     
    2727   PUBLIC   sbc_oce_alloc   ! routine called in sbcmod.F90 
    2828   PUBLIC   sbc_tau2wnd     ! routine called in several sbc modules 
    29     
     29 
    3030   !!---------------------------------------------------------------------- 
    3131   !!           Namelist for the Ocean Surface Boundary Condition 
     
    4545   LOGICAL , PUBLIC ::   ln_dm2dc       !: Daily mean to Diurnal Cycle short wave (qsr) 
    4646   LOGICAL , PUBLIC ::   ln_rnf         !: runoffs / runoff mouths 
    47    LOGICAL , PUBLIC ::   ln_ssr         !: Sea Surface restoring on SST and/or SSS       
     47   LOGICAL , PUBLIC ::   ln_ssr         !: Sea Surface restoring on SST and/or SSS 
    4848   LOGICAL , PUBLIC ::   ln_apr_dyn     !: Atmospheric pressure forcing used on dynamics (ocean & ice) 
    4949   INTEGER , PUBLIC ::   nn_ice         !: flag for ice in the surface boundary condition (=0/1/2/3) 
     
    5151   !                                             !: =F levitating ice (no presure effect) with mass and salt exchanges 
    5252   !                                             !: =T embedded sea-ice (pressure effect + mass and salt exchanges) 
    53    INTEGER , PUBLIC ::   nn_components  !: flag for sbc module (including sea-ice) coupling mode (see component definition below)  
    54    INTEGER , PUBLIC ::   nn_fwb         !: FreshWater Budget:  
    55    !                                             !:  = 0 unchecked  
     53   INTEGER , PUBLIC ::   nn_components  !: flag for sbc module (including sea-ice) coupling mode (see component definition below) 
     54   INTEGER , PUBLIC ::   nn_fwb         !: FreshWater Budget: 
     55   !                                             !:  = 0 unchecked 
    5656   !                                             !:  = 1 global mean of e-p-r set to zero at each nn_fsbc time step 
    5757   !                                             !:  = 2 annual global mean of e-p-r set to zero 
     
    8181   INTEGER , PUBLIC, PARAMETER ::   jp_purecpl = 5        !: Pure ocean-atmosphere Coupled formulation 
    8282   INTEGER , PUBLIC, PARAMETER ::   jp_none    = 6        !: for OPA when doing coupling via SAS module 
    83     
    84    !!---------------------------------------------------------------------- 
    85    !!           Stokes drift parametrization definition  
     83 
     84   !!---------------------------------------------------------------------- 
     85   !!           Stokes drift parametrization definition 
    8686   !!---------------------------------------------------------------------- 
    8787   INTEGER , PUBLIC, PARAMETER ::   jp_breivik_2014 = 0     !: Breivik  2014: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 
    88    INTEGER , PUBLIC, PARAMETER ::   jp_li_2017      = 1     !: Li et al 2017: Stokes drift based on Phillips spectrum (Breivik 2016)  
    89                                                             !  with depth averaged profile 
    90    INTEGER , PUBLIC, PARAMETER ::   jp_peakfr       = 2     !: Li et al 2017: using the peak wave number read from wave model instead  
    91                                                             !  of the inverse depth scale 
     88   INTEGER , PUBLIC, PARAMETER ::   jp_li_2017      = 1     !: Li et al 2017: Stokes drift based on Phillips spectrum (Breivik 2016) 
     89   !  with depth averaged profile 
     90   INTEGER , PUBLIC, PARAMETER ::   jp_peakfr       = 2     !: Li et al 2017: using the peak wave number read from wave model instead 
     91   !  of the inverse depth scale 
    9292   LOGICAL , PUBLIC            ::   ll_st_bv2014  = .FALSE. !  logical indicator, .true. if Breivik 2014 parameterisation is active. 
    9393   LOGICAL , PUBLIC            ::   ll_st_li2017  = .FALSE. !  logical indicator, .true. if Li 2017 parameterisation is active. 
     
    9898   !!           component definition 
    9999   !!---------------------------------------------------------------------- 
    100    INTEGER , PUBLIC, PARAMETER ::   jp_iam_nemo = 0      !: Initial single executable configuration  
    101                                                          !  (no internal OASIS coupling) 
     100   INTEGER , PUBLIC, PARAMETER ::   jp_iam_nemo = 0      !: Initial single executable configuration 
     101   !  (no internal OASIS coupling) 
    102102   INTEGER , PUBLIC, PARAMETER ::   jp_iam_opa  = 1      !: Multi executable configuration - OPA component 
    103                                                          !  (internal OASIS coupling) 
     103   !  (internal OASIS coupling) 
    104104   INTEGER , PUBLIC, PARAMETER ::   jp_iam_sas  = 2      !: Multi executable configuration - SAS component 
    105                                                          !  (internal OASIS coupling) 
     105   !  (internal OASIS coupling) 
    106106   !!---------------------------------------------------------------------- 
    107107   !!              Ocean Surface Boundary Condition fields 
     
    112112   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2] 
    113113   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2] 
    114    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum              !: module of sea surface stress (at T-point)    [N/m2]  
     114   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum              !: module of sea surface stress (at T-point)    [N/m2] 
    115115   !! wndm is used compute surface gases exchanges in ice-free ocean or leads 
    116116   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s] 
    117    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rhoa              !: air density at "rn_zu" m above the sea       [kg/m3] !LB 
     117   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rhoa              !: air density at "rn_zu" m above the sea       [kg/m3] 
    118118   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsr               !: sea heat flux:     solar                     [W/m2] 
    119119   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qns    , qns_b    !: sea heat flux: non solar                     [W/m2] 
     
    124124   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp_tot           !: total E-P over ocean and ice                 [Kg/m2/s] 
    125125   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmflx            !: freshwater budget: freezing/melting          [Kg/m2/s] 
    126    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff                                 [Kg/m2/s]   
    127    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwficb , fwficb_b !: iceberg melting                              [Kg/m2/s]   
     126   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff                                 [Kg/m2/s] 
     127   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwficb , fwficb_b !: iceberg melting                              [Kg/m2/s] 
    128128   !! 
    129129   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  sbc_tsc, sbc_tsc_b  !: sbc content trend                      [K.m/s] jpi,jpj,jpts 
     
    138138 
    139139   !!--------------------------------------------------------------------- 
    140    !! ABL Vertical Domain size   
     140   !! ABL Vertical Domain size 
    141141   !!--------------------------------------------------------------------- 
    142142   INTEGER , PUBLIC            ::   jpka   = 2     !: ABL number of vertical levels (default definition) 
     
    154154   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sss_m     !: mean (nn_fsbc time-step) surface sea salinity            [psu] 
    155155   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssh_m     !: mean (nn_fsbc time-step) sea surface height                [m] 
     156   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk_m     !: mean (nn_fsbc time-step) SKIN surface sea temp.      [Celsius] 
    156157   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_m     !: mean (nn_fsbc time-step) sea surface layer thickness       [m] 
    157158   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frq_m     !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-] 
     
    175176      ! 
    176177      ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , taum(jpi,jpj) ,     & 
    177          &      vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) )  
    178          ! 
     178         &      vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) ) 
     179      ! 
    179180      ALLOCATE( qns_tot(jpi,jpj) , qns  (jpi,jpj) , qns_b(jpi,jpj),        & 
    180181         &      qsr_tot(jpi,jpj) , qsr  (jpi,jpj) ,                        & 
    181182         &      emp    (jpi,jpj) , emp_b(jpi,jpj) ,                        & 
    182183         &      sfx    (jpi,jpj) , sfx_b(jpi,jpj) , emp_tot(jpi,jpj), fmmflx(jpi,jpj), STAT=ierr(2) ) 
    183          ! 
     184      ! 
    184185      ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,  & 
    185186         &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) ,  & 
    186187         &      fwficb  (jpi,jpj), fwficb_b(jpi,jpj), STAT=ierr(3) ) 
    187          ! 
     188      ! 
    188189      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     & 
    189          &      atm_co2(jpi,jpj) ,                                        & 
     190         &      atm_co2(jpi,jpj) , tsk_m(jpi,jpj) ,                       & 
    190191         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      & 
    191192         &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 
     
    203204      !!--------------------------------------------------------------------- 
    204205      !!                    ***  ROUTINE sbc_tau2wnd  *** 
    205       !!                    
    206       !! ** Purpose : Estimation of wind speed as a function of wind stress    
     206      !! 
     207      !! ** Purpose : Estimation of wind speed as a function of wind stress 
    207208      !! 
    208209      !! ** Method  : |tau|=rhoa*Cd*|U|^2 
     
    215216      INTEGER  ::   ji, jj                ! dummy indices 
    216217      !!--------------------------------------------------------------------- 
    217       zcoef = 0.5 / ( zrhoa * zcdrag )  
     218      zcoef = 0.5 / ( zrhoa * zcdrag ) 
    218219      DO jj = 2, jpjm1 
    219220         DO ji = fs_2, fs_jpim1   ! vect. opt. 
    220             ztx = utau(ji-1,jj  ) + utau(ji,jj)  
    221             zty = vtau(ji  ,jj-1) + vtau(ji,jj)  
     221            ztx = utau(ji-1,jj  ) + utau(ji,jj) 
     222            zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
    222223            ztau = SQRT( ztx * ztx + zty * zty ) 
    223224            wndm(ji,jj) = SQRT ( ztau * zcoef ) * tmask(ji,jj,1) 
  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/sbcblk.F90

    r12193 r12199  
    178178      ! 
    179179      !                             !** read bulk namelist 
     180      REWIND( numnam_ref )                !* Namelist namsbc_blk in reference namelist : bulk parameters 
    180181      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) 
    181182901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' ) 
    182183      ! 
     184      REWIND( numnam_cfg )                !* Namelist namsbc_blk in configuration namelist : bulk parameters 
    183185      READ  ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 ) 
    184186902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' ) 
     
    399401      ! Sanity/consistence test on humidity at first time step to detect potential screw-up: 
    400402      IF( kt == nit000 ) THEN 
    401          WRITE(numout,*) '' 
     403         IF(lwp) WRITE(numout,*) '' 
    402404#if defined key_agrif 
    403          WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ===' 
     405         IF(lwp) WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ===' 
    404406#else 
    405407         ztmp = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain 
     
    415417            END SELECT 
    416418            IF(ztmp < 0._wp) THEN 
    417                WRITE(numout,'("   Mean humidity value found on proc #",i5.5," is: ",f)') narea, ztmp 
     419               IF (lwp) WRITE(numout,'("   Mean humidity value found on proc #",i5.5," is: ",f)') narea, ztmp 
    418420               CALL ctl_stop( 'STOP', 'Something is wrong with air humidity!!!', & 
    419421                  &   ' ==> check the unit in your input files'       , & 
     
    422424            END IF 
    423425         END IF 
    424          WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 
     426         IF(lwp) WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 
    425427#endif 
    426          WRITE(numout,*) '' 
     428         IF(lwp) WRITE(numout,*) '' 
    427429      END IF !IF( kt == nit000 ) 
    428430      !                                            ! compute the surface ocean fluxes using bulk formulea 
     
    432434            &                sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m,       &   !   <<= in 
    433435            &                sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1),   &   !   <<= in (wl/cs) 
    434             &                zssq, zcd_du, zsen, zevp )                              !   =>> out 
     436            &                tsk_m, zssq, zcd_du, zsen, zevp )                       !   =>> out 
    435437 
    436438         CALL blk_oce_2(     sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),   &   !   <<= in 
    437439            &                sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1),   &   !   <<= in 
    438             &                sf(jp_snow)%fnow(:,:,1), sst_m,                     &   !   <<= in 
     440            &                sf(jp_snow)%fnow(:,:,1), tsk_m,                     &   !   <<= in 
    439441            &                zsen, zevp )                                            !   <=> in out 
    440442      ENDIF 
     
    470472 
    471473   SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, &  ! inp 
    472       &              pslp , pst   , pu   , pv,    &  ! inp 
    473       &              pqsr , pqlw  ,               &  ! inp 
    474       &              pssq , pcd_du, psen , pevp   )  ! out 
     474      &                  pslp , pst   , pu   , pv,        &  ! inp 
     475      &                  pqsr , pqlw  ,                   &  ! inp 
     476      &                  ptsk, pssq , pcd_du, psen , pevp   )  ! out 
    475477      !!--------------------------------------------------------------------- 
    476478      !!                     ***  ROUTINE blk_oce_1  *** 
     
    494496      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   ptair  ! potential temperature at T-points        [Kelvin] 
    495497      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pslp   ! sea-level pressure                       [Pa] 
    496       REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pst    ! surface temperature                      [Celcius] 
     498      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pst    ! surface temperature                      [Celsius] 
    497499      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pu     ! surface current at U-point (i-component) [m/s] 
    498500      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pv     ! surface current at V-point (j-component) [m/s] 
    499501      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqsr   ! 
    500502      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqlw   ! 
     503      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   ptsk   ! skin temp. (or SST if CS & WL not used)  [Celsius] 
    501504      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pssq   ! specific humidity at pst                 [kg/kg] 
    502505      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pcd_du ! Cd x |dU| at T-points                    [m/s] 
     
    507510      REAL(wp) ::   zztmp                ! local variable 
    508511      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
    509       REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
    510512      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    511513      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
     
    519521      ! 
    520522      ! local scalars ( place there for vector optimisation purposes) 
    521       zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
     523      !                           ! Temporary conversion from Celcius to Kelvin (and set minimum value far above 0 K) 
     524      ptsk(:,:) = pst(:,:) + rt0  ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!) 
    522525 
    523526      ! ----------------------------------------------------------------------------- ! 
     
    566569 
    567570      ! specific humidity at SST 
    568       pssq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), pslp(:,:) ) 
     571      pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) ) 
    569572 
    570573      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    571          zztmp1(:,:) = zst(:,:) 
     574         !! Backup "bulk SST" and associated spec. hum. 
     575         zztmp1(:,:) = ptsk(:,:) 
    572576         zztmp2(:,:) = pssq(:,:) 
    573577      ENDIF 
     
    608612 
    609613      CASE( np_NCAR      ) 
    610          CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm,                              & 
     614         CALL turb_ncar    ( rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm,                              & 
    611615            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    612616 
    613617      CASE( np_COARE_3p0 ) 
    614          CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     618         CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
    615619            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
    616620            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
    617621 
    618622      CASE( np_COARE_3p6 ) 
    619          CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     623         CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
    620624            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
    621625            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
    622626 
    623627      CASE( np_ECMWF     ) 
    624          CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
     628         CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
    625629            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
    626630            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
     
    632636 
    633637      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    634          !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and pssq: 
    635          WHERE ( fr_i < 0.001_wp ) 
    636             ! zst and pssq have been updated by cool-skin/warm-layer scheme and we keep it!!! 
    637             zst(:,:)  =  zst(:,:)*tmask(:,:,1) 
    638             pssq(:,:) = pssq(:,:)*tmask(:,:,1) 
    639          ELSEWHERE 
    640             ! we forget about the update... 
    641             zst(:,:)  = zztmp1(:,:) !#LB: using what we backed up before skin-algo 
    642             pssq(:,:) = zztmp2(:,:) !#LB:  "   "   " 
     638         !! ptsk and pssq have been updated!!! 
     639         !! 
     640         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of ptsk and pssq: 
     641         WHERE ( fr_i(:,:) > 0.001_wp ) 
     642            ! sea-ice present, we forget about the update, using what we backed up before call to turb_*() 
     643            ptsk(:,:) = zztmp1(:,:) 
     644            pssq(:,:) = zztmp2(:,:) 
    643645         END WHERE 
    644646      END IF 
     
    669671         END DO 
    670672      ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation 
    671          CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 
     673         CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 
    672674            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         & 
    673675            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                 & 
     
    707709         ENDIF 
    708710         ! 
    709       ENDIF 
    710       ! 
     711      ENDIF !IF( ln_abl ) 
     712       
     713      ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1)  ! Back to Celsius 
     714             
     715      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     716         CALL iom_put( "t_skin" ,  ptsk        )  ! T_skin in Celsius 
     717         CALL iom_put( "dt_skin" , ptsk - pst  )  ! T_skin - SST temperature difference... 
     718      ENDIF 
     719 
    711720      IF(ln_ctl) THEN 
    712721         CALL prt_ctl( tab2d_1=pevp  , clinfo1=' blk_oce_1: pevp   : ' ) 
     
    719728 
    720729   SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec,   &   ! <<= in 
    721       &          psnow, pst , psen, pevp     )   ! <<= in 
     730      &                  psnow, ptsk, psen, pevp     )   ! <<= in 
    722731      !!--------------------------------------------------------------------- 
    723732      !!                     ***  ROUTINE blk_oce_2  *** 
     
    740749      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pprec 
    741750      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psnow 
    742       REAL(wp), INTENT(in), DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius] 
     751      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptsk   ! SKIN surface temperature   [Celsius] 
    743752      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psen 
    744753      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pevp 
     
    746755      INTEGER  ::   ji, jj               ! dummy loop indices 
    747756      REAL(wp) ::   zztmp,zz1,zz2,zz3    ! local variable 
    748       REAL(wp), DIMENSION(jpi,jpj) ::   zqlw              ! long wave and sensible heat fluxes 
     757      REAL(wp), DIMENSION(jpi,jpj) ::   ztskk             ! skin temp. in Kelvin  
     758      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw              ! long wave and sensible heat fluxes       
    749759      REAL(wp), DIMENSION(jpi,jpj) ::   zqla              ! latent heat fluxes and evaporation 
    750       REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
    751760      !!--------------------------------------------------------------------- 
    752761      ! 
    753762      ! local scalars ( place there for vector optimisation purposes) 
    754       zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    755  
    756  
     763 
     764      ztskk(:,:) = ptsk(:,:) + rt0  ! => ptsk in Kelvin rather than Celsius 
     765       
    757766      ! ----------------------------------------------------------------------------- ! 
    758767      !     III    Net longwave radiative FLUX                                        ! 
     
    760769 
    761770      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 
    762       !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 
    763       zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux 
    764  
    765       !  Turbulent fluxes over ocean 
    766       ! ----------------------------- 
     771      !! (ztskk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 
     772      zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*ztskk(:,:)*ztskk(:,:)*ztskk(:,:)*ztskk(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux 
     773 
     774      !  Latent flux over ocean 
     775      ! ----------------------- 
    767776 
    768777      ! use scalar version of L_vap() for AGRIF compatibility 
    769778      DO jj = 1, jpj 
    770779         DO ji = 1, jpi 
    771             zqla(ji,jj) = -1._wp * L_vap( zst(ji,jj) ) * pevp(ji,jj)    ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 
     780            zqla(ji,jj) = - L_vap( ztskk(ji,jj) ) * pevp(ji,jj)    ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 
    772781         ENDDO 
    773782      ENDDO 
     
    788797      qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:)                   &   ! Downward Non Solar 
    789798         &     - psnow(:,:) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip 
    790          &     - pevp(:,:) * pst(:,:) * rcp                          &   ! remove evap heat content at SST !LB??? pst is Celsius !? 
     799         &     - pevp(:,:) * ptsk(:,:) * rcp                         &   ! remove evap heat content at SST 
    791800         &     + ( pprec(:,:) - psnow(:,:) ) * rn_pfac               &   ! add liquid precip heat content at Tair 
    792801         &     * ( ptair(:,:) - rt0 ) * rcp                          & 
     
    815824         CALL iom_put( "qsr_oce"  ,   qsr  )               ! output downward solar heat over the ocean 
    816825         CALL iom_put( "qt_oce"   ,   qns+qsr )            ! output total downward heat over the ocean 
    817       ENDIF 
    818       ! 
    819       IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    820          CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius 
    821          CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference... 
    822826      ENDIF 
    823827      ! 
     
    11051109 
    11061110      IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN 
    1107          ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )  
     1111         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) 
    11081112         IF( iom_use('evap_ao_cea'  ) )  CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average) 
    11091113         IF( iom_use('hflx_evap_cea') )  CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) )   ! heat flux from evap (cell average) 
     
    11141118      ENDIF 
    11151119      IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN 
    1116           WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) ;   ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 
    1117           ELSEWHERE                             ;   ztmp(:,:) = rcp * sst_m(:,:)     
    1118           ENDWHERE 
    1119           ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus )  
    1120           IF( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average) 
    1121           IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean) 
    1122           IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice) 
     1120         WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) 
     1121            ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 
     1122         ELSEWHERE 
     1123            ztmp(:,:) = rcp * sst_m(:,:) 
     1124         ENDWHERE 
     1125         ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus ) 
     1126         IF( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average) 
     1127         IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean) 
     1128         IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice) 
    11231129      ENDIF 
    11241130      ! 
  • NEMO/branches/2019/dev_r11943_MERGE_2019/tests/STATION_ASF/EXPREF/context_nemo.xml

    r11833 r12199  
    11<!-- 
    2  ==============================================================================================  
     2    ============================================================================================== 
    33    NEMO context 
    4 ==============================================================================================  
     4    ============================================================================================== 
    55--> 
    66<context id="nemo"> 
    7 <!-- $id$ --> 
    8 <!-- Fields definition --> 
    9     <field_definition src="./field_def_nemo-oce.xml"/>   <!--  NEMO ocean dynamics                     --> 
     7  <variable_definition> 
     8    <!-- Year/Month/Day of time origin for NetCDF files; defaults to 1800-01-01 --> 
     9    <variable id="ref_year"  type="int"> 1970 </variable> 
     10    <variable id="ref_month" type="int"> 01 </variable> 
     11    <variable id="ref_day"   type="int"> 01 </variable> 
     12  </variable_definition> 
     13  <!-- --> 
     14  <!-- Fields definition --> 
     15  <field_definition src="./field_def_nemo-oce.xml"/>   <!--  NEMO ocean dynamics                     --> 
    1016 
    11 <!-- Files definition --> 
    12     <file_definition src="./file_def_nemo-oce.xml"/>     <!--  NEMO ocean dynamics                     --> 
    13     <!--  
    14 ============================================================================================================ 
    15 = grid definition = = DO NOT CHANGE = 
    16 ============================================================================================================ 
    17     --> 
    18      
    19     <axis_definition> 
    20       <axis id="deptht" long_name="Vertical T levels" unit="m" positive="down" /> 
    21       <axis id="depthu" long_name="Vertical U levels" unit="m" positive="down" /> 
    22       <axis id="depthv" long_name="Vertical V levels" unit="m" positive="down" /> 
    23       <axis id="depthw" long_name="Vertical W levels" unit="m" positive="down" /> 
    24       <axis id="nfloat" long_name="Float number"      unit="-"                 /> 
    25       <axis id="icbcla"  long_name="Iceberg class"      unit="1"               /> 
    26       <axis id="ncatice" long_name="Ice category"       unit="1"               /> 
    27       <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"            /> 
    28       <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"            /> 
    29     </axis_definition> 
    30   
    31     <domain_definition src="./domain_def_nemo.xml"/> 
    32    
    33     <grid_definition src="./grid_def_nemo.xml"/> 
    34      
     17  <!-- Files definition --> 
     18  <file_definition src="./file_def_nemo-oce.xml"/>     <!--  NEMO ocean dynamics                     --> 
     19  <!-- 
     20      ============================================================================================================ 
     21      = grid definition = = DO NOT CHANGE = 
     22      ============================================================================================================ 
     23  --> 
     24 
     25  <axis_definition> 
     26    <axis id="deptht" long_name="Vertical T levels" unit="m" positive="down" /> 
     27    <axis id="depthu" long_name="Vertical U levels" unit="m" positive="down" /> 
     28    <axis id="depthv" long_name="Vertical V levels" unit="m" positive="down" /> 
     29    <axis id="depthw" long_name="Vertical W levels" unit="m" positive="down" /> 
     30    <axis id="nfloat" long_name="Float number"      unit="-"                 /> 
     31    <axis id="icbcla"  long_name="Iceberg class"      unit="1"               /> 
     32    <axis id="ncatice" long_name="Ice category"       unit="1"               /> 
     33    <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"            /> 
     34    <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"            /> 
     35  </axis_definition> 
     36 
     37  <domain_definition src="./domain_def_nemo.xml"/> 
     38 
     39  <grid_definition src="./grid_def_nemo.xml"/> 
     40 
    3541</context> 
  • NEMO/branches/2019/dev_r11943_MERGE_2019/tests/STATION_ASF/EXPREF/namelist_coare3p6-noskin_cfg

    r12126 r12199  
    4343      cn_ocerst_out   = 'restart_oce'   !  suffix of ocean restart name (output) 
    4444      cn_ocerst_outdir = './'         !  directory in which to write output ocean restarts 
    45    ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model 
    4645   nn_istate   =       0   !  output the initial state (1) or not (0) 
    4746   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 
  • NEMO/branches/2019/dev_r11943_MERGE_2019/tests/STATION_ASF/EXPREF/namelist_coare3p6_cfg

    r12126 r12199  
    4343      cn_ocerst_out   = 'restart_oce'   !  suffix of ocean restart name (output) 
    4444      cn_ocerst_outdir = './'         !  directory in which to write output ocean restarts 
    45    ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model 
    4645   nn_istate   =       0   !  output the initial state (1) or not (0) 
    4746   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 
  • NEMO/branches/2019/dev_r11943_MERGE_2019/tests/STATION_ASF/EXPREF/namelist_ecmwf-noskin_cfg

    r12126 r12199  
    4343      cn_ocerst_out   = 'restart_oce'   !  suffix of ocean restart name (output) 
    4444      cn_ocerst_outdir = './'         !  directory in which to write output ocean restarts 
    45    ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model 
    4645   nn_istate   =       0   !  output the initial state (1) or not (0) 
    4746   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 
  • NEMO/branches/2019/dev_r11943_MERGE_2019/tests/STATION_ASF/EXPREF/namelist_ecmwf_cfg

    r12126 r12199  
    4343      cn_ocerst_out   = 'restart_oce'   !  suffix of ocean restart name (output) 
    4444      cn_ocerst_outdir = './'         !  directory in which to write output ocean restarts 
    45    ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model 
    4645   nn_istate   =       0   !  output the initial state (1) or not (0) 
    4746   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 
  • NEMO/branches/2019/dev_r11943_MERGE_2019/tests/STATION_ASF/EXPREF/namelist_ncar_cfg

    r12126 r12199  
    4343      cn_ocerst_out   = 'restart_oce'   !  suffix of ocean restart name (output) 
    4444      cn_ocerst_outdir = './'         !  directory in which to write output ocean restarts 
    45    ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model 
    4645   nn_istate   =       0   !  output the initial state (1) or not (0) 
    4746   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 
  • NEMO/branches/2019/dev_r11943_MERGE_2019/tests/STATION_ASF/MY_SRC/usrdef_nam.F90

    r11930 r12199  
    6060      !!---------------------------------------------------------------------- 
    6161      ! 
    62       REWIND( numnam_cfg )          ! Namelist namusr_def (exist in namelist_cfg only) 
    6362      READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 ) 
    6463902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist' ) 
Note: See TracChangeset for help on using the changeset viewer.