- Timestamp:
- 2019-12-11T12:09:17+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/ABL/sbcabl.F90
r12154 r12179 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 76 & ln_smth_pblh 77 77 !!--------------------------------------------------------------------- 78 78 … … 87 87 IF(lwm) WRITE( numond, namsbc_abl ) 88 88 ! 89 ! Check ABL mixing length option 89 ! Check ABL mixing length option 90 90 IF( nn_amxl < 0 .OR. nn_amxl > 2 ) & 91 & CALL ctl_stop( 'abl_init : bad flag, nn_amxl must be 0, 1 or 2 ' ) 92 ! 93 ! 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 94 94 IF( nn_dyn_restore < 0 .OR. nn_dyn_restore > 2 ) & 95 & 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 ' ) 96 96 97 97 !!--------------------------------------------------------------------- 98 98 !! Control prints 99 !!--------------------------------------------------------------------- 99 !!--------------------------------------------------------------------- 100 100 IF(lwp) THEN ! Control print (other namelist variable) 101 101 WRITE(numout,*) … … 106 106 IF(ln_geos_winds) THEN 107 107 ln_geos_winds = .FALSE. 108 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.)' 109 109 END IF 110 110 ELSE IF( ln_geos_winds ) THEN 111 WRITE(numout,*) ' ABL -- winds forced by geostrophic winds' 111 WRITE(numout,*) ' ABL -- winds forced by geostrophic winds' 112 112 ELSE 113 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' 114 114 END IF 115 115 ! 116 116 SELECT CASE ( nn_dyn_restore ) 117 CASE ( 0 ) 118 WRITE(numout,*) ' ABL -- No restoring for ABL winds' 117 CASE ( 0 ) 118 WRITE(numout,*) ' ABL -- No restoring for ABL winds' 119 119 CASE ( 1 ) 120 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 ' 121 121 CASE ( 2 ) 122 WRITE(numout,*) ' ABL -- Restoring of ABL winds activated everywhere ' 122 WRITE(numout,*) ' ABL -- Restoring of ABL winds activated everywhere ' 123 123 END SELECT 124 124 ! 125 IF( ln_smth_pblh ) WRITE(numout,*) ' ABL -- Smoothing of PBL height is activated' 126 125 IF( ln_smth_pblh ) WRITE(numout,*) ' ABL -- Smoothing of PBL height is activated' 126 ! 127 127 ENDIF 128 128 129 129 !!--------------------------------------------------------------------- 130 130 !! Convert nudging coefficient from hours to 1/sec 131 !!--------------------------------------------------------------------- 131 !!--------------------------------------------------------------------- 132 132 zcff = 1._wp / 3600._wp 133 133 rn_ldyn_min = zcff / rn_ldyn_min 134 rn_ldyn_max = zcff / rn_ldyn_max 134 rn_ldyn_max = zcff / rn_ldyn_max 135 135 rn_ltra_min = zcff / rn_ltra_min 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 ) 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 ) 142 142 id = iom_varid( inum, 'e3t_abl', kdimsz=idimsz, kndims=indims, lduld=lluldl ) 143 jpka = idimsz(indims - COUNT( (/lluldl/) ) ) 143 jpka = idimsz(indims - COUNT( (/lluldl/) ) ) 144 144 jpkam1 = jpka - 1 145 145 … … 152 152 153 153 #if ! defined key_iomput 154 154 IF( dia_wri_alloc_abl() /= 0 ) CALL ctl_stop( 'STOP', 'abl_init : unable to allocate arrays' ) 155 155 #endif 156 156 … … 158 158 WRITE(numout,*) 159 159 WRITE(numout,*) ' sbc_abl_init : ABL Reference vertical grid' 160 WRITE(numout,*) ' ~~~~~~~' 160 WRITE(numout,*) ' ~~~~~~~' 161 161 WRITE(numout, "(9x,' level ght_abl ghw_abl e3t_abl e3w_abl ')" ) 162 162 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, ght_abl(jk), ghw_abl(jk), e3t_abl(jk), e3w_abl(jk), jk = 1, jpka ) … … 164 164 165 165 !!--------------------------------------------------------------------- 166 !! Check TKE closure parameters 167 !!--------------------------------------------------------------------- 166 !! Check TKE closure parameters 167 !!--------------------------------------------------------------------- 168 168 rn_Sch = rn_ce / rn_cm 169 169 mxl_min = (avm_bak / rn_cm) / sqrt( tke_min ) … … 172 172 WRITE(numout,*) 173 173 WRITE(numout,*) ' abl_zdf_tke : ABL TKE turbulent closure' 174 WRITE(numout,*) ' ~~~~~~~~~~~' 174 WRITE(numout,*) ' ~~~~~~~~~~~' 175 175 IF(nn_amxl==0) WRITE(numout,*) 'Deardorff 80 length-scale ' 176 176 IF(nn_amxl==1) WRITE(numout,*) 'length-scale based on the distance to the PBL height ' 177 177 WRITE(numout,*) ' Minimum value of atmospheric TKE = ',tke_min,' m^2 s^-2' 178 178 WRITE(numout,*) ' Minimum value of atmospheric mixing length = ',mxl_min,' m' 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 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 183 183 END IF 184 184 185 185 !!------------------------------------------------------------------------------------------- 186 186 !! Compute parameters to build the vertical profile for the nudging term (used in abl_stp()) 187 !!------------------------------------------------------------------------------------------- 187 !!------------------------------------------------------------------------------------------- 188 188 zcff1 = 1._wp / ( jp_bmax - jp_bmin )**3 189 189 ! for active tracers … … 192 192 jp_alp1_tra = -6._wp * zcff1 * jp_bmax * jp_bmin * ( rn_ltra_max - rn_ltra_min ) 193 193 jp_alp0_tra = zcff1 * ( rn_ltra_max * jp_bmin*jp_bmin * (3._wp*jp_bmax - jp_bmin) & 194 & - 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) ) 195 195 ! for dynamics 196 196 jp_alp3_dyn = -2._wp * zcff1 * ( rn_ldyn_max - rn_ldyn_min ) … … 198 198 jp_alp1_dyn = -6._wp * zcff1 * jp_bmax * jp_bmin * ( rn_ldyn_max - rn_ldyn_min ) 199 199 jp_alp0_dyn = zcff1 * ( rn_ldyn_max * jp_bmin*jp_bmin * (3._wp*jp_bmax - jp_bmin) & 200 & - 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) ) 201 201 202 202 jp_pblh_min = ghw_abl( 4) / jp_bmin !<-- at least 3 grid points at the bottom have value rn_ltra_min … … 229 229 & + jp_alp1_tra * jp_bmin + jp_alp0_tra ) * rdt_abl 230 230 zcff1 = ( jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2 & 231 & + jp_alp1_tra * jp_bmax + jp_alp0_tra ) * rdt_abl 231 & + jp_alp1_tra * jp_bmax + jp_alp0_tra ) * rdt_abl 232 232 IF(lwp) THEN 233 233 WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff … … 244 244 !!------------------------------------------------------------------------------------------- 245 245 !! Initialize Coriolis frequency, equatorial restoring and land/sea mask 246 !!------------------------------------------------------------------------------------------- 246 !!------------------------------------------------------------------------------------------- 247 247 fft_abl(:,:) = 2._wp * omega * SIN( rad * gphit(:,:) ) 248 248 249 249 ! Equatorial restoring 250 250 IF( nn_dyn_restore == 1 ) THEN 251 zcff = 2._wp * omega * SIN( rad * 90._wp ) !++ fmax 251 zcff = 2._wp * omega * SIN( rad * 90._wp ) !++ fmax 252 252 rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff) / zcff ) )**8 253 253 !!GS: alternative shape … … 260 260 msk_abl(:,:) = tmask(:,:,1) 261 261 262 !!------------------------------------------------------------------------------------------- 262 !!------------------------------------------------------------------------------------------- 263 263 264 264 ! initialize 2D bulk fields AND 3D abl data … … 275 275 CALL fld_read( nit000, nn_fsbc, sf ) ! input fields provided at the first time-step 276 276 277 278 277 u_abl(:,:,:,nt_n ) = sf(jp_wndi)%fnow(:,:,:) 278 v_abl(:,:,:,nt_n ) = sf(jp_wndj)%fnow(:,:,:) 279 279 tq_abl(:,:,:,nt_n,jp_ta) = sf(jp_tair)%fnow(:,:,:) 280 280 tq_abl(:,:,:,nt_n,jp_qa) = sf(jp_humi)%fnow(:,:,:) 281 281 282 282 tke_abl(:,:,:,nt_n ) = tke_min 283 283 avm_abl(:,:,: ) = avm_bak 284 284 avt_abl(:,:,: ) = avt_bak 285 mxl_abl(:,:,: ) = mxl_min 286 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 287 287 u_abl (:,:,:,nt_a ) = 0._wp 288 288 v_abl (:,:,:,nt_a ) = 0._wp … … 292 292 293 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) 294 294 295 295 END SUBROUTINE sbc_abl_init 296 296 … … 301 301 !! 302 302 !! ** Purpose : provide the momentum, heat and freshwater fluxes at 303 !! the ocean surface from an ABL calculation at each oceanic time step 304 !! 305 !! ** Method : 303 !! the ocean surface from an ABL calculation at each oceanic time step 304 !! 305 !! ** Method : 306 306 !! - Pre-compute part of turbulent fluxes in blk_oce_1 307 !! - Perform 1 time-step of the ABL model 307 !! - Perform 1 time-step of the ABL model 308 308 !! - Finalize flux computation in blk_oce_2 309 309 !! … … 322 322 #if defined key_si3 323 323 REAL(wp), DIMENSION(jpi,jpj) :: zssqi, zcd_dui, zseni, zevpi 324 #endif 324 #endif 325 325 INTEGER :: jbak, jbak_dta, ji, jj 326 326 !!--------------------------------------------------------------------- … … 329 329 !! 1 - Read Atmospheric 3D data for large-scale forcing 330 330 !!------------------------------------------------------------------------------------------- 331 331 332 332 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 333 333 334 334 !!------------------------------------------------------------------------------------------- 335 335 !! 2 - Compute Cd x ||U||, Ch x ||U||, Ce x ||U||, and SSQ using now fields 336 !!------------------------------------------------------------------------------------------- 336 !!------------------------------------------------------------------------------------------- 337 337 338 338 CALL blk_oce_1( kt, u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in … … 340 340 & sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m , & ! <<= in 341 341 & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) , & ! <<= in 342 & zssq, zcd_du, zsen, zevp )! =>> out343 342 & tsk_m, zssq, zcd_du, zsen, zevp ) ! =>> out 343 344 344 #if defined key_si3 345 345 CALL blk_ice_1( u_abl(:,:,2,nt_n ), v_abl(:,:,2,nt_n ), & ! <<= in … … 347 347 & sf(jp_slp)%fnow(:,:,1) , u_ice, v_ice, tm_su , & ! <<= in 348 348 & pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui ) ! <<= out 349 #endif 349 #endif 350 350 351 351 !!------------------------------------------------------------------------------------------- 352 352 !! 3 - Advance ABL variables from now (n) to after (n+1) 353 !!------------------------------------------------------------------------------------------- 354 355 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 356 356 & sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:), & ! <<= in 357 357 & sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:), & ! <<= in … … 360 360 & zcd_du, zsen, zevp, & ! <=> in/out 361 361 & wndm, utau, vtau, taum & ! =>> out 362 #if defined key_si3 362 #if defined key_si3 363 363 & , tm_su, u_ice, v_ice, zssqi, zcd_dui & ! <<= in 364 364 & , zseni, zevpi, wndm_ice, ato_i & ! <<= in 365 365 & , utau_ice, vtau_ice & ! =>> out 366 #endif 366 #endif 367 367 & ) 368 368 !!------------------------------------------------------------------------------------------- 369 !! 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 370 370 !! time swap is done in abl_stp 371 !!------------------------------------------------------------------------------------------- 371 !!------------------------------------------------------------------------------------------- 372 372 373 373 CALL blk_oce_2( tq_abl(:,:,2,nt_n,jp_ta), & 374 374 & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1), & 375 375 & sf(jp_prec)%fnow(:,:,1) , sf(jp_snow)%fnow(:,:,1), & 376 & sst_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 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 380 380 381 381 #if defined key_si3 382 382 ! Avoid a USE abl in icesbc module 383 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) 384 #endif 384 #endif 385 385 386 386 END SUBROUTINE sbc_abl
Note: See TracChangeset
for help on using the changeset viewer.