- Timestamp:
- 2019-07-29T16:47:01+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.
