Changeset 11363
- Timestamp:
- 2019-07-29T16:47:01+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D
- Files:
-
- 8 edited
-
cfgs/SHARED/field_def_nemo-ice.xml (modified) (1 diff)
-
cfgs/SHARED/field_def_nemo-oce.xml (modified) (5 diffs)
-
cfgs/SHARED/grid_def_nemo.xml (modified) (1 diff)
-
cfgs/SHARED/namelist_ref (modified) (2 diffs)
-
src/ABL/ablmod.F90 (modified) (19 diffs)
-
src/ABL/par_abl.F90 (modified) (2 diffs)
-
src/ABL/sbcabl.F90 (modified) (3 diffs)
-
src/OCE/DIA/diawri.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-ice.xml
r10911 r11363 90 90 <field id="qemp_ice" long_name="Downward Heat Flux from E-P over ice" unit="W/m2" /> 91 91 <field id="albedo" long_name="Mean albedo over sea ice and ocean" unit="" /> 92 <field id="Cd_ice" long_name="Momentum turbulent exchange coefficient" unit="" /> 93 <field id="Ch_ice" long_name="Heat turbulent exchange coefficient" unit="" /> 92 94 93 95 <!-- trends --> -
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-oce.xml
r11306 r11363 24 24 25 25 <field id="toce" long_name="temperature" standard_name="sea_water_potential_temperature" unit="degC" grid_ref="grid_T_3D"/> 26 <field id="toce_e3t" long_name="temperature (thickness weighted)" unit="degC" grid_ref="grid_T_3D" > toce * e3t </field 26 <field id="toce_e3t" long_name="temperature (thickness weighted)" unit="degC" grid_ref="grid_T_3D" > toce * e3t </field> 27 27 <field id="soce" long_name="salinity" standard_name="sea_water_practical_salinity" unit="1e-3" grid_ref="grid_T_3D"/> 28 <field id="soce_e3t" long_name="salinity (thickness weighted)" unit="1e-3" grid_ref="grid_T_3D" > soce * e3t </field 28 <field id="soce_e3t" long_name="salinity (thickness weighted)" unit="1e-3" grid_ref="grid_T_3D" > soce * e3t </field> 29 29 30 30 <!-- t-eddy viscosity coefficients (ldfdyn) --> … … 33 33 34 34 <field id="sst" long_name="sea surface temperature" standard_name="sea_surface_temperature" unit="degC" /> 35 <field id="sst2" long_name="square of sea surface temperature" standard_name="square_of_sea_surface_temperature" unit="degC2" > sst * sst </field 35 <field id="sst2" long_name="square of sea surface temperature" standard_name="square_of_sea_surface_temperature" unit="degC2" > sst * sst </field> 36 36 <field id="sstmax" long_name="max of sea surface temperature" field_ref="sst" operation="maximum" /> 37 37 <field id="sstmin" long_name="min of sea surface temperature" field_ref="sst" operation="minimum" /> … … 45 45 46 46 <field id="sss" long_name="sea surface salinity" standard_name="sea_surface_salinity" unit="1e-3" /> 47 <field id="sss2" long_name="square of sea surface salinity" unit="1e-6" > sss * sss </field 47 <field id="sss2" long_name="square of sea surface salinity" unit="1e-6" > sss * sss </field> 48 48 <field id="sssmax" long_name="max of sea surface salinity" field_ref="sss" operation="maximum" /> 49 49 <field id="sssmin" long_name="min of sea surface salinity" field_ref="sss" operation="minimum" /> … … 54 54 55 55 <field id="ssh" long_name="sea surface height" standard_name="sea_surface_height_above_geoid" unit="m" /> 56 <field id="ssh2" long_name="square of sea surface height" standard_name="square_of_sea_surface_height_above_geoid" unit="m2" > ssh * ssh </field 56 <field id="ssh2" long_name="square of sea surface height" standard_name="square_of_sea_surface_height_above_geoid" unit="m2" > ssh * ssh </field> 57 57 <field id="wetdep" long_name="wet depth" standard_name="wet_depth" unit="m" /> 58 58 <field id="sshmax" long_name="max of sea surface height" field_ref="ssh" operation="maximum" /> … … 352 352 <field_group id="ABL" > <!-- time step automaticaly defined based on nn_fsbc --> 353 353 <!-- variables available with ABL on atmospheric T grid--> 354 <field_group id="grid_ABL " grid_ref="grid_TA_3D" >354 <field_group id="grid_ABL3D" grid_ref="grid_TA_3D" > 355 355 <field id="u_abl" long_name="ABL i-horizontal velocity" standard_name="abl_x_velocity" unit="m/s" /> 356 356 <field id="v_abl" long_name="ABL j-horizontal velocity" standard_name="abl_y_velocity" unit="m/s" /> 357 <field id="t_abl" long_name="ABL potential temperature" standard_name="abl_theta" unit="K" /> 358 <field id="q_abl" long_name="ABL specific humidity" standard_name="abl_qspe" unit="kg/kg" /> 359 <field id="pblh" long_name="ABL height" standard_name="abl_height" unit="m" grid_ref="grid_T_2D" /> 357 <field id="t_abl" long_name="ABL potential temperature" standard_name="abl_theta" unit="K" /> 358 <field id="q_abl" long_name="ABL specific humidity" standard_name="abl_qspe" unit="kg/kg" /> 359 <field id="coeft" long_name="ABL nudging coefficient" standard_name="coeft" unit="" /> 360 </field_group> 361 <field_group id="grid_ABL2D" grid_ref="grid_TA_2D" > 362 <field id="pblh" long_name="ABL height" standard_name="abl_height" unit="m" /> 360 363 </field_group> 361 364 </field_group> <!-- ABL --> -
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/grid_def_nemo.xml
r11275 r11363 59 59 60 60 <!-- ABL grid definition --> 61 <grid id="grid_TA_2D"> 62 <domain domain_ref="grid_T" /> 63 </grid> 61 64 <grid id="grid_TA_3D"> 62 65 <domain domain_ref="grid_T" /> 63 66 <axis id="ght_abl" /> 64 </grid>65 <grid id="grid_TA_2D">66 <domain domain_ref="grid_T" />67 67 </grid> 68 68 <grid id="grid_WA_3D"> -
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ref
r11334 r11363 283 283 &namsbc_abl ! Atmospheric Boundary Layer formulation (ln_abl = T) 284 284 !----------------------------------------------------------------------- 285 cn_dir = ' atm_forcing/' ! root directory for the location of the ABL grid file285 cn_dir = './' ! root directory for the location of the ABL grid file 286 286 cn_dom = 'dom_cfg_abl.nc' 287 287 ln_hpgls_frc = .false. … … 301 301 rn_Rod = 0.15 ! c0 in RMCA17 mixing length formulation (not yet implemented) 302 302 rn_Ric = 0.139 ! Critical Richardson number (to compute PBL height and diffusivities) 303 ln_smth_pblh = . true.! Smoothing of PBL height with a 2d Hanning filter303 ln_smth_pblh = .false. ! Smoothing of PBL height with a 2d Hanning filter 304 304 / 305 305 !----------------------------------------------------------------------- -
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90
r11348 r11363 153 153 DO jk = 3, jpkam1 154 154 DO ji = 1, jpi ! vector opt. 155 z_elem_a( ji, jk ) = - rdt * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal156 z_elem_c( ji, jk ) = - rdt * Avt_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal155 z_elem_a( ji, jk ) = - rdt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal 156 z_elem_c( ji, jk ) = - rdt_abl * Avt_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 157 157 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal 158 158 END DO … … 162 162 ! Neumann at the bottom 163 163 z_elem_a( ji, 2 ) = 0._wp 164 z_elem_c( ji, 2 ) = - rdt * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 )164 z_elem_c( ji, 2 ) = - rdt_abl * Avt_abl( ji, jj, 2 ) / e3w_abl( 2 ) 165 165 ! Homogeneous Neumann at the top 166 z_elem_a( ji, jpka ) = - rdt * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )166 z_elem_a( ji, jpka ) = - rdt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka ) 167 167 z_elem_c( ji, jpka ) = 0._wp 168 168 z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) … … 185 185 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj) 186 186 #endif 187 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1188 tq_abl ( ji, jj, 2, nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rdt * zztmp2187 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1 188 tq_abl ( ji, jj, 2, nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rdt_abl * zztmp2 189 189 tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra ) 190 190 END DO … … 197 197 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj) 198 198 #endif 199 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1200 tq_abl ( ji, jj, 2, nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rdt * zztmp2199 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1 200 tq_abl ( ji, jj, 2, nt_a, jtra ) = e3t_abl( 2 ) * tq_abl ( ji, jj, 2, nt_n, jtra ) + rdt_abl * zztmp2 201 201 tq_abl ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl ( ji, jj, jpka, nt_n, jtra ) 202 202 END DO … … 244 244 DO jj = 1, jpj 245 245 DO ji = 1, jpi 246 zcff = ( fft_abl(ji,jj) * rdt )*( fft_abl(ji,jj) * rdt) ! (f dt)**2246 zcff = ( fft_abl(ji,jj) * rdt_abl )*( fft_abl(ji,jj) * rdt_abl ) ! (f dt)**2 247 247 248 248 u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 249 249 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*u_abl( ji, jj, jk, nt_n ) & 250 & + rdt * fft_abl(ji, jj) * v_abl ( ji , jj , jk, nt_n ) ) &250 & + rdt_abl * fft_abl(ji, jj) * v_abl ( ji , jj , jk, nt_n ) ) & 251 251 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 252 252 253 253 v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( & 254 254 & (1._wp-gamma_Cor*(1._wp-gamma_Cor)*zcff)*v_abl( ji, jj, jk, nt_n ) & 255 & - rdt * fft_abl(ji, jj) * u_abl ( ji , jj, jk, nt_n ) ) &255 & - rdt_abl * fft_abl(ji, jj) * u_abl ( ji , jj, jk, nt_n ) ) & 256 256 & / (1._wp + gamma_Cor*gamma_Cor*zcff) 257 257 END DO … … 267 267 DO ji = 1, jpi 268 268 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) & 269 & - rdt * e3t_abl(jk) * fft_abl(ji , jj) * pgv_dta(ji ,jj ,jk)269 & - rdt_abl * e3t_abl(jk) * fft_abl(ji , jj) * pgv_dta(ji ,jj ,jk) 270 270 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) & 271 & + rdt * e3t_abl(jk) * fft_abl(ji, jj ) * pgu_dta(ji ,jj ,jk)271 & + rdt_abl * e3t_abl(jk) * fft_abl(ji, jj ) * pgu_dta(ji ,jj ,jk) 272 272 END DO 273 273 END DO … … 280 280 DO jk = 1, jpka 281 281 DO ji = 1, jpi 282 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgu_dta(ji,jj,jk)283 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rdt * e3t_abl(jk) * pgv_dta(ji,jj,jk)282 u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rdt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk) 283 v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rdt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk) 284 284 ENDDO 285 285 ENDDO … … 298 298 DO jk = 3, jpkam1 299 299 DO ji = 1, jpi 300 z_elem_a( ji, jk ) = - rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal301 z_elem_c( ji, jk ) = - rdt * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal300 z_elem_a( ji, jk ) = - rdt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal 301 z_elem_c( ji, jk ) = - rdt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 302 302 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal 303 303 END DO … … 307 307 !++ Surface boundary condition 308 308 z_elem_a( ji, 2 ) = 0._wp 309 z_elem_c( ji, 2 ) = - rdt * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )309 z_elem_c( ji, 2 ) = - rdt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 310 310 ! 311 311 zztmp1 = pcd_du(ji, jj) … … 316 316 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 317 317 #endif 318 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1319 u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rdt * zztmp2318 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1 319 u_abl( ji, jj, 2, nt_a ) = u_abl( ji, jj, 2, nt_a ) + rdt_abl * zztmp2 320 320 321 321 !++ Top Neumann B.C. 322 !z_elem_a( ji, jpka ) = - 0.5_wp * rdt * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka )322 !z_elem_a( ji, jpka ) = - 0.5_wp * rdt_abl * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka ) 323 323 !z_elem_c( ji, jpka ) = 0._wp 324 324 !z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) … … 365 365 DO jk = 3, jpkam1 366 366 DO ji = 1, jpi 367 z_elem_a( ji, jk ) = -rdt * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal368 z_elem_c( ji, jk ) = -rdt * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal367 z_elem_a( ji, jk ) = -rdt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 ) ! lower-diagonal 368 z_elem_c( ji, jk ) = -rdt_abl * Avm_abl( ji, jj, jk ) / e3w_abl( jk ) ! upper-diagonal 369 369 z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) ! diagonal 370 370 END DO … … 374 374 !++ Surface boundary condition 375 375 z_elem_a( ji, 2 ) = 0._wp 376 z_elem_c( ji, 2 ) = - rdt * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )376 z_elem_c( ji, 2 ) = - rdt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 ) 377 377 ! 378 378 zztmp1 = pcd_du(ji, jj) … … 383 383 zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 384 384 #endif 385 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1386 v_abl( ji, jj, 2, nt_a ) = v_abl( ji, jj, 2, nt_a ) + rdt * zztmp2385 z_elem_b( ji, 2 ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt_abl * zztmp1 386 v_abl( ji, jj, 2, nt_a ) = v_abl( ji, jj, 2, nt_a ) + rdt_abl * zztmp2 387 387 !++ Top Neumann B.C. 388 !z_elem_a( ji, jpka ) = -rdt * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )388 !z_elem_a( ji, jpka ) = -rdt_abl * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka ) 389 389 !z_elem_c( ji, jpka ) = 0._wp 390 390 !z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka ) … … 434 434 DO jj = 2, jpj 435 435 DO ji = 2, jpi 436 zcff1 437 zsig 438 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) )436 zcff1 = pblh( ji, jj ) 437 zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) ) 438 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 439 439 zmsk = msk_abl(ji,jj) 440 440 zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2 & 441 441 & + jp_alp1_dyn * zsig + jp_alp0_dyn 442 zcff = (1._wp-zmsk) + zmsk * rdt * zcff2! zcff = 1 for masked points443 444 442 zcff = (1._wp-zmsk) + zmsk * zcff2 * rdt ! zcff = 1 for masked points 443 ! rdt = rdt_abl / nn_fsbc 444 zcff = zcff * rest_eq(ji,jj) ; z_cft( ji, jj, jk ) = zcff 445 445 446 446 u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * u_abl( ji, jj, jk, nt_a ) & … … 448 448 v_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) * v_abl( ji, jj, jk, nt_a ) & 449 449 & + zcff * pv_dta( ji, jj, jk ) 450 450 END DO 451 451 END DO 452 452 !------------- … … 460 460 DO jj = 1,jpj 461 461 DO ji = 1,jpi 462 zcff1 = pblh( ji, jj ) 463 zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) ) 464 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 465 zmsk = msk_abl(ji,jj) 466 zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2 & 467 & + jp_alp1_tra * zsig + jp_alp0_tra 468 zcff = (1._wp-zmsk) + zmsk * rdt * zcff2 ! zcff = 1 for masked points 462 zcff1 = pblh( ji, jj ) 463 zsig = ght_abl(jk) / MAX( jp_pblh_min, MIN( jp_pblh_max, zcff1 ) ) 464 zsig = MIN( jp_bmax , MAX( zsig, jp_bmin) ) 465 zmsk = msk_abl(ji,jj) 466 zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2 & 467 & + jp_alp1_tra * zsig + jp_alp0_tra 468 zcff = (1._wp-zmsk) + zmsk * zcff2 * rdt ! zcff = 1 for masked points 469 ! rdt = rdt_abl / nn_fsbc 469 470 tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta ) & 470 471 & + zcff * pt_dta( ji, jj, jk ) … … 483 484 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 484 485 ! 485 CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a), 'T', -1., v_abl(:,:,:,nt_a), 'T', -1., kfillmode = jpfillnothing )486 CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1., tq_abl(:,:,:,nt_a,jp_qa), 'T',1., kfillmode = jpfillnothing ) ! ++++ this should not be needed...486 CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a ), 'T', -1., v_abl(:,:,:,nt_a ), 'T', -1., kfillmode = jpfillnothing ) 487 CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1., tq_abl(:,:,:,nt_a,jp_qa), 'T', 1., kfillmode = jpfillnothing ) ! ++++ this should not be needed... 487 488 ! 488 489 CALL iom_put ( 'u_abl', u_abl(:,:,2:jpka,nt_a )) … … 490 491 CALL iom_put ( 't_abl', tq_abl(:,:,2:jpka,nt_a,jp_ta)) 491 492 CALL iom_put ( 'q_abl', tq_abl(:,:,2:jpka,nt_a,jp_qa)) 492 !CALL iom_put ( 'coeft', z_cft(:,:,2:jpka ))493 CALL iom_put ( 'coeft', z_cft(:,:,2:jpka )) 493 494 CALL iom_put ( 'pblh', pblh(:,: )) 494 495 ! … … 676 677 zbuoy = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk ) 677 678 678 z_elem_a( ji, jk ) = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk ) ! lower-diagonal679 z_elem_c( ji, jk ) = - 0.5_wp * rdt * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal679 z_elem_a( ji, jk ) = - 0.5_wp * rdt_abl * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk-1 ) ) / e3t_abl( jk ) ! lower-diagonal 680 z_elem_c( ji, jk ) = - 0.5_wp * rdt_abl * rn_Sch * ( Avm_abl( ji, jj, jk )+Avm_abl( ji, jj, jk+1 ) ) / e3t_abl( jk+1 ) ! upper-diagonal 680 681 IF( (zbuoy + zshear) .gt. 0.) THEN ! Patankar trick to avoid negative values of TKE 681 682 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 682 & + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) ! diagonal683 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt * ( zbuoy + zshear ) ) ! right-hand-side683 & + e3w_abl(jk) * rdt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) ! diagonal 684 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt_abl * ( zbuoy + zshear ) ) ! right-hand-side 684 685 ELSE 685 686 z_elem_b( ji, jk ) = e3w_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk ) & 686 & + e3w_abl(jk) * rdt * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) & ! diagonal687 & - e3w_abl(jk) * rdt * zbuoy688 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt * zshear ) ! right-hand-side687 & + e3w_abl(jk) * rdt_abl * rn_Ceps * sqrt(tke_abl( ji, jj, jk, nt_n )) / mxl_abl(ji,jj,jk) & ! diagonal 688 & - e3w_abl(jk) * rdt_abl * zbuoy 689 tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rdt_abl * zshear ) ! right-hand-side 689 690 END IF 690 691 END DO -
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/par_abl.F90
r11322 r11363 19 19 INTEGER , PUBLIC, PARAMETER :: jp_qa = 2 !: indice for humidity 20 20 INTEGER , PUBLIC, PARAMETER :: jptime = 2 !: number of time indices stored in memory 21 21 22 !!--------------------------------------------------------------------- 22 23 !! namABL Namelist options … … 69 70 ! parameter for the semi-implicit treatment of Coriolis term 70 71 REAL(wp), PUBLIC, PARAMETER :: gamma_Cor = 0.55_wp 72 ! ABL timestep 73 REAL(wp), PUBLIC :: rdt_abl 71 74 72 75 !!--------------------------------------------------------------------- -
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90
r11360 r11363 200 200 jp_pblh_max = ghw_abl(jpka-3) / jp_bmax !<-- at least 3 grid points at the top have value rn_ltra_max 201 201 202 ! ABL timestep 203 rdt_abl = nn_fsbc * rdt 204 202 205 ! Check parameters for dynamics 203 zcff = jp_alp3_dyn * jp_bmin**3 + jp_alp2_dyn * jp_bmin**2 &204 & + jp_alp1_dyn * jp_bmin + jp_alp0_dyn205 zcff1 = jp_alp3_dyn * jp_bmax**3 + jp_alp2_dyn * jp_bmax**2 &206 & + jp_alp1_dyn * jp_bmax + jp_alp0_dyn206 zcff = ( jp_alp3_dyn * jp_bmin**3 + jp_alp2_dyn * jp_bmin**2 & 207 & + jp_alp1_dyn * jp_bmin + jp_alp0_dyn ) * rdt_abl 208 zcff1 = ( jp_alp3_dyn * jp_bmax**3 + jp_alp2_dyn * jp_bmax**2 & 209 & + jp_alp1_dyn * jp_bmax + jp_alp0_dyn ) * rdt_abl 207 210 IF(lwp) THEN 208 211 IF(nn_dyn_restore > 0) THEN 209 WRITE(numout,*) ' ABL Minimum value for dynamics restoring = ',zcff * rdt210 WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1 * rdt212 WRITE(numout,*) ' ABL Minimum value for dynamics restoring = ',zcff 213 WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1 211 214 ! Check that restoring coefficients are between 0 and 1 212 IF( zcff1 * rdt > 1._wp .OR. zcff1 * rdt < 0._wp ) & 215 !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp ) & 216 IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp ) & 213 217 & CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_max' ) 214 IF( zcff * rdt > 1._wp .OR. zcff * rdt < 0._wp ) & 218 !IF( zcff > 1._wp .OR. zcff < 0._wp ) & 219 IF( zcff > nn_fsbc .OR. zcff < 0._wp ) & 215 220 & CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_min' ) 216 IF( zcff * rdt > zcff1 * rdt )&221 IF( zcff > zcff1 ) & 217 222 & CALL ctl_stop( 'abl_init : rn_ldyn_max should be smaller than rn_ldyn_min' ) 218 223 END IF … … 220 225 221 226 ! Check parameters for active tracers 222 zcff = jp_alp3_tra * jp_bmin**3 + jp_alp2_tra * jp_bmin**2 &223 & + jp_alp1_tra * jp_bmin + jp_alp0_tra224 zcff1 = jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2 &225 & + jp_alp1_tra * jp_bmax + jp_alp0_tra227 zcff = ( jp_alp3_tra * jp_bmin**3 + jp_alp2_tra * jp_bmin**2 & 228 & + jp_alp1_tra * jp_bmin + jp_alp0_tra ) * rdt_abl 229 zcff1 = ( jp_alp3_tra * jp_bmax**3 + jp_alp2_tra * jp_bmax**2 & 230 & + jp_alp1_tra * jp_bmax + jp_alp0_tra ) * rdt_abl 226 231 IF(lwp) THEN 227 WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff * rdt228 WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1 * rdt232 WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff 233 WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1 229 234 ! Check that restoring coefficients are between 0 and 1 230 IF( zcff1 * rdt > 1._wp .OR. zcff1 * rdt < 0._wp ) & 235 !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp ) & 236 IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp ) & 231 237 & CALL ctl_stop( 'abl_init : wrong value for rn_ltra_max' ) 232 IF( zcff * rdt > 1._wp .OR. zcff * rdt < 0._wp ) & 238 !IF( zcff > 1._wp .OR. zcff < 0._wp ) & 239 IF( zcff > nn_fsbc .OR. zcff < 0._wp ) & 233 240 & CALL ctl_stop( 'abl_init : wrong value for rn_ltra_min' ) 234 IF( zcff * rdt > zcff1 * rdt )&241 IF( zcff > zcff1 ) & 235 242 & CALL ctl_stop( 'abl_init : rn_ltra_max should be smaller than rn_ltra_min' ) 236 243 END IF … … 244 251 IF( nn_dyn_restore == 1 ) THEN 245 252 zcff = 2._wp * omega * SIN( rad * 90._wp ) !++ fmax 246 rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff)/zcff ) )**8 253 rest_eq(:,:) = SIN( 0.5_wp*rpi*( (fft_abl(:,:) - zcff) / zcff ) )**8 254 !!GS: alternative shape 255 !rest_eq(:,:) = SIN( 0.5_wp*rpi*(zcff - ABS(ff_t(:,:))) / (zcff - 3.e-5) )**8 256 !WHERE(ABS(ff_t(:,:)).LE.3.e-5) rest_eq(:,:) = 1._wp 247 257 ELSE 248 258 rest_eq(:,:) = 1._wp -
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90
r11360 r11363 654 654 CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm 655 655 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 656 !656 ! 657 657 IF( ln_abl ) THEN 658 ! !!! nid_A : 3D 659 CALL histdef( nid_A, "t_abl", "Potential Temperature" , "K" , & ! t_abl 658 CALL histdef( nid_A, "t_abl", "Potential Temperature" , "K" , & ! t_abl 660 659 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 661 660 CALL histdef( nid_A, "q_abl", "Humidity" , "kg/kg" , & ! q_abl … … 677 676 & jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 678 677 #endif 679 CALL histend( nid_A, snc4chunks=snc4set ) 680 ! 681 ENDIF 682 ! 678 CALL histend( nid_A, snc4chunks=snc4set ) 679 ENDIF 680 ! 683 681 IF( ln_icebergs ) THEN 684 682 CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , & … … 838 836 CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction 839 837 CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed 840 !838 ! 841 839 IF( ln_abl ) THEN 842 843 844 845 846 END DO 847 840 ALLOCATE( zw3d_abl(jpi,jpj,jpka) ) 841 IF( ln_mskland ) THEN 842 DO jk=1,jpka 843 zw3d_abl(:,:,jk) = tmask(:,:,1) 844 END DO 845 ELSE 848 846 zw3d_abl(:,:,:) = 1._wp 849 847 ENDIF 850 851 852 853 854 855 856 857 848 CALL histwrite( nid_A, "pblh" , it, pblh(:,:) *zw3d_abl(:,:,1 ), ndim_hA, ndex_hA ) ! pblh 849 CALL histwrite( nid_A, "u_abl" , it, u_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! u_abl 850 CALL histwrite( nid_A, "v_abl" , it, v_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! v_abl 851 CALL histwrite( nid_A, "t_abl" , it, tq_abl (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! t_abl 852 CALL histwrite( nid_A, "q_abl" , it, tq_abl (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! q_abl 853 CALL histwrite( nid_A, "tke_abl", it, tke_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! tke_abl 854 CALL histwrite( nid_A, "avm_abl", it, avm_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avm_abl 855 CALL histwrite( nid_A, "avt_abl", it, avt_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avt_abl 858 856 #if defined key_si3 859 857 CALL histwrite( nid_A, "oce_frac" , it, ato_i(:,:) , ndim_hA, ndex_hA ) ! ato_i 860 858 #endif 861 862 863 !859 DEALLOCATE(zw3d_abl) 860 ENDIF 861 ! 864 862 IF( ln_icebergs ) THEN 865 863 !
Note: See TracChangeset
for help on using the changeset viewer.
