Changeset 12199 for NEMO/branches/2019/dev_r11943_MERGE_2019/src/ABL
- Timestamp:
- 2019-12-12T09:30:08+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/ABL/sbcabl.F90
r12193 r12199 2 2 !!====================================================================== 3 3 !! *** 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 6 6 !!===================================================================== 7 7 !! History : 4.0 ! 2019-03 (F. Lemarié & G. Samson) Original code … … 9 9 10 10 !!---------------------------------------------------------------------- 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 13 13 !! fluxes over ocean via the ABL model 14 14 !!---------------------------------------------------------------------- … … 35 35 USE ice , ONLY : u_ice, v_ice, tm_su, ato_i ! ato_i = total open water fractional area 36 36 USE sbc_ice, ONLY : wndm_ice, utau_ice, vtau_ice 37 #endif 37 #endif 38 38 #if ! defined key_iomput 39 39 USE diawri , ONLY : dia_wri_alloc_abl … … 59 59 !! 60 60 !! ** 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 63 63 !! 64 64 !!---------------------------------------------------------------------- … … 74 74 & rn_ldyn_min , rn_ldyn_max, rn_ltra_min, rn_ltra_max, & 75 75 & 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 79 80 READ ( numnam_ref, namsbc_abl, IOSTAT = ios, ERR = 901 ) 80 81 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in reference namelist' ) 81 82 ! 83 REWIND( numnam_cfg ) ! Namelist namsbc_abl in configuration namelist : ABL parameters 82 84 READ ( numnam_cfg, namsbc_abl, IOSTAT = ios, ERR = 902 ) 83 85 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_abl in configuration namelist' ) … … 85 87 IF(lwm) WRITE( numond, namsbc_abl ) 86 88 ! 87 ! Check ABL mixing length option 89 ! Check ABL mixing length option 88 90 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 92 94 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 ' ) 94 96 95 97 !!--------------------------------------------------------------------- 96 98 !! Control prints 97 !!--------------------------------------------------------------------- 99 !!--------------------------------------------------------------------- 98 100 IF(lwp) THEN ! Control print (other namelist variable) 99 101 WRITE(numout,*) … … 104 106 IF(ln_geos_winds) THEN 105 107 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.)' 107 109 END IF 108 110 ELSE IF( ln_geos_winds ) THEN 109 WRITE(numout,*) ' ABL -- winds forced by geostrophic winds' 111 WRITE(numout,*) ' ABL -- winds forced by geostrophic winds' 110 112 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' 112 114 END IF 113 115 ! 114 116 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' 117 119 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 ' 119 121 CASE ( 2 ) 120 WRITE(numout,*) ' ABL -- Restoring of ABL winds activated everywhere ' 122 WRITE(numout,*) ' ABL -- Restoring of ABL winds activated everywhere ' 121 123 END SELECT 122 124 ! 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 ! 125 127 ENDIF 126 128 127 129 !!--------------------------------------------------------------------- 128 130 !! Convert nudging coefficient from hours to 1/sec 129 !!--------------------------------------------------------------------- 131 !!--------------------------------------------------------------------- 130 132 zcff = 1._wp / 3600._wp 131 133 rn_ldyn_min = zcff / rn_ldyn_min 132 rn_ldyn_max = zcff / rn_ldyn_max 134 rn_ldyn_max = zcff / rn_ldyn_max 133 135 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 ) 140 142 id = iom_varid( inum, 'e3t_abl', kdimsz=idimsz, kndims=indims, lduld=lluldl ) 141 jpka = idimsz(indims - COUNT( (/lluldl/) ) ) 143 jpka = idimsz(indims - COUNT( (/lluldl/) ) ) 142 144 jpkam1 = jpka - 1 143 145 … … 150 152 151 153 #if ! defined key_iomput 152 154 IF( dia_wri_alloc_abl() /= 0 ) CALL ctl_stop( 'STOP', 'abl_init : unable to allocate arrays' ) 153 155 #endif 154 156 … … 156 158 WRITE(numout,*) 157 159 WRITE(numout,*) ' sbc_abl_init : ABL Reference vertical grid' 158 WRITE(numout,*) ' ~~~~~~~' 160 WRITE(numout,*) ' ~~~~~~~' 159 161 WRITE(numout, "(9x,' level ght_abl ghw_abl e3t_abl e3w_abl ')" ) 160 162 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, ght_abl(jk), ghw_abl(jk), e3t_abl(jk), e3w_abl(jk), jk = 1, jpka ) … … 162 164 163 165 !!--------------------------------------------------------------------- 164 !! Check TKE closure parameters 165 !!--------------------------------------------------------------------- 166 !! Check TKE closure parameters 167 !!--------------------------------------------------------------------- 166 168 rn_Sch = rn_ce / rn_cm 167 169 mxl_min = (avm_bak / rn_cm) / sqrt( tke_min ) … … 170 172 WRITE(numout,*) 171 173 WRITE(numout,*) ' abl_zdf_tke : ABL TKE turbulent closure' 172 WRITE(numout,*) ' ~~~~~~~~~~~' 174 WRITE(numout,*) ' ~~~~~~~~~~~' 173 175 IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale ' 174 176 IF(nn_amxl==1) WRITE(numout,*) 'length-scale based on the distance to the PBL height ' 175 177 WRITE(numout,*) ' Minimum value of atmospheric TKE = ',tke_min,' m^2 s^-2' 176 178 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 181 183 END IF 182 184 183 185 !!------------------------------------------------------------------------------------------- 184 186 !! Compute parameters to build the vertical profile for the nudging term (used in abl_stp()) 185 !!------------------------------------------------------------------------------------------- 187 !!------------------------------------------------------------------------------------------- 186 188 zcff1 = 1._wp / ( jp_bmax - jp_bmin )**3 187 189 ! for active tracers … … 190 192 jp_alp1_tra = -6._wp * zcff1 * jp_bmax * jp_bmin * ( rn_ltra_max - rn_ltra_min ) 191 193 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) ) 193 195 ! for dynamics 194 196 jp_alp3_dyn = -2._wp * zcff1 * ( rn_ldyn_max - rn_ldyn_min ) … … 196 198 jp_alp1_dyn = -6._wp * zcff1 * jp_bmax * jp_bmin * ( rn_ldyn_max - rn_ldyn_min ) 197 199 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) ) 199 201 200 202 jp_pblh_min = ghw_abl( 4) / jp_bmin !<-- at least 3 grid points at the bottom have value rn_ltra_min … … 227 229 & + jp_alp1_tra * jp_bmin + jp_alp0_tra ) * rdt_abl 228 230 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 230 232 IF(lwp) THEN 231 233 WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff … … 242 244 !!------------------------------------------------------------------------------------------- 243 245 !! Initialize Coriolis frequency, equatorial restoring and land/sea mask 244 !!------------------------------------------------------------------------------------------- 246 !!------------------------------------------------------------------------------------------- 245 247 fft_abl(:,:) = 2._wp * omega * SIN( rad * gphit(:,:) ) 246 248 247 249 ! Equatorial restoring 248 250 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 250 252 rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff) / zcff ) )**8 251 253 !!GS: alternative shape … … 258 260 msk_abl(:,:) = tmask(:,:,1) 259 261 260 !!------------------------------------------------------------------------------------------- 262 !!------------------------------------------------------------------------------------------- 261 263 262 264 ! initialize 2D bulk fields AND 3D abl data … … 273 275 CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step 274 276 275 276 277 u_abl(:,:,:,nt_n ) = sf(jp_wndi)%fnow(:,:,:) 278 v_abl(:,:,:,nt_n ) = sf(jp_wndj)%fnow(:,:,:) 277 279 tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:) 278 280 tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:) 279 281 280 282 tke_abl(:,:,:,nt_n ) = tke_min 281 283 avm_abl(:,:,: ) = avm_bak 282 284 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 285 287 u_abl (:,:,:,nt_a ) = 0._wp 286 288 v_abl (:,:,:,nt_a ) = 0._wp … … 290 292 291 293 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 293 295 END SUBROUTINE sbc_abl_init 294 296 … … 299 301 !! 300 302 !! ** 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 : 304 306 !! - 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 306 308 !! - Finalize flux computation in blk_oce_2 307 309 !! … … 320 322 #if defined key_si3 321 323 REAL(wp), DIMENSION(jpi,jpj) :: zssqi, zcd_dui, zseni, zevpi 322 #endif 324 #endif 323 325 INTEGER :: jbak, jbak_dta, ji, jj 324 326 !!--------------------------------------------------------------------- … … 327 329 !! 1 - Read Atmospheric 3D data for large-scale forcing 328 330 !!------------------------------------------------------------------------------------------- 329 331 330 332 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 331 333 332 334 !!------------------------------------------------------------------------------------------- 333 335 !! 2 - Compute Cd x ||U||, Ch x ||U||, Ce x ||U||, and SSQ using now fields 334 !!------------------------------------------------------------------------------------------- 336 !!------------------------------------------------------------------------------------------- 335 337 336 338 CALL blk_oce_1( kt, u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in … … 338 340 & sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m , & ! <<= in 339 341 & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) , & ! <<= in 340 & zssq, zcd_du, zsen, zevp )! =>> out341 342 & tsk_m, zssq, zcd_du, zsen, zevp ) ! =>> out 343 342 344 #if defined key_si3 343 345 CALL blk_ice_1( u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in … … 345 347 & sf(jp_slp)%fnow(:,:,1) , u_ice, v_ice, tm_su , & ! <<= in 346 348 & pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui ) ! <<= out 347 #endif 349 #endif 348 350 349 351 !!------------------------------------------------------------------------------------------- 350 352 !! 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, & ! <<= in353 !!------------------------------------------------------------------------------------------- 354 355 CALL abl_stp( kt, tsk_m, ssu_m, ssv_m, zssq, & ! <<= in 354 356 & sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:), & ! <<= in 355 357 & sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:), & ! <<= in … … 358 360 & zcd_du, zsen, zevp, & ! <=> in/out 359 361 & wndm, utau, vtau, taum & ! =>> out 360 #if defined key_si3 362 #if defined key_si3 361 363 & , tm_su, u_ice, v_ice, zssqi, zcd_dui & ! <<= in 362 364 & , zseni, zevpi, wndm_ice, ato_i & ! <<= in 363 365 & , utau_ice, vtau_ice & ! =>> out 364 #endif 366 #endif 365 367 & ) 366 368 !!------------------------------------------------------------------------------------------- 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 368 370 !! time swap is done in abl_stp 369 !!------------------------------------------------------------------------------------------- 371 !!------------------------------------------------------------------------------------------- 370 372 371 373 CALL blk_oce_2( tq_abl(:,:,2,nt_n,jp_ta), & 372 374 & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1), & 373 375 & 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 378 380 379 381 #if defined key_si3 380 382 ! Avoid a USE abl in icesbc module 381 383 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 383 385 384 386 END SUBROUTINE sbc_abl
Note: See TracChangeset
for help on using the changeset viewer.
