- Timestamp:
- 2021-06-15T16:39:31+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/ticket2669_isf_fluxes/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/ticket2669_isf_fluxes/src
- Property svn:mergeinfo deleted
-
NEMO/branches/2021/ticket2669_isf_fluxes/src/OCE/ZDF/zdfgls.F90
r14834 r14994 26 26 USE zdfmxl ! mixed layer 27 27 USE sbcwave , ONLY : hsw ! significant wave height 28 #if defined key_si3 29 USE ice, ONLY: hm_i, h_i 30 #endif 31 #if defined key_cice 32 USE sbc_ice, ONLY: h_i 33 #endif 28 34 ! 29 35 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 51 57 LOGICAL :: ln_length_lim ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 52 58 LOGICAL :: ln_sigpsi ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing 59 INTEGER :: nn_mxlice ! type of scaling under sea-ice (=0/1/2/3) 53 60 INTEGER :: nn_bc_surf ! surface boundary condition (=0/1) 54 61 INTEGER :: nn_bc_bot ! bottom boundary condition (=0/1) … … 208 215 zhsro(:,:) = rn_hsro 209 216 CASE ( 1 ) ! Standard Charnock formula 210 zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(A2D(nn_hls)) , rn_hsro ) 217 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 218 zhsro(ji,jj) = MAX( rsbc_zs1 * ustar2_surf(ji,jj) , rn_hsro ) 219 END_2D 211 220 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 212 221 !!gm faster coding : the 2 comment lines should be used … … 222 231 ! 223 232 ! adapt roughness where there is sea ice 224 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 225 zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * rn_hsri )*tmask(ji,jj,1) + & 226 & (1._wp - tmask(ji,jj,1))*rn_hsro 227 END_2D 233 SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice 234 ! 235 CASE( 1 ) ! scaling with constant sea-ice roughness 236 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 237 zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * rn_hsri )*tmask(ji,jj,1) + (1._wp - tmask(ji,jj,1))*rn_hsro 238 END_2D 239 ! 240 CASE( 2 ) ! scaling with mean sea-ice thickness 241 #if defined key_si3 242 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 243 zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * hm_i(ji,jj) )*tmask(ji,jj,1) + (1._wp - tmask(ji,jj,1))*rn_hsro 244 END_2D 245 #endif 246 ! 247 CASE( 3 ) ! scaling with max sea-ice thickness 248 #if defined key_si3 || defined key_cice 249 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 250 zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * MAXVAL(h_i(ji,jj,:)) )*tmask(ji,jj,1) + (1._wp - tmask(ji,jj,1))*rn_hsro 251 END_2D 252 #endif 253 ! 254 END SELECT 228 255 ! 229 256 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !== Compute dissipation rate ==! … … 327 354 END_2D 328 355 ! 356 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 357 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 358 IF( mikt(ji,jj) > 1 )THEN 359 itop = mikt(ji,jj) ! k top w-point 360 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 361 ! ! mask at the 362 ! ocean surface 363 ! points 364 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 365 ! 366 ! Dirichlet condition applied at: 367 ! top level (itop) & Just below it (itopp1) 368 zd_lw(ji,jj,itop) = 0._wp ; zd_lw(ji,jj,itopp1) = 0._wp 369 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 370 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = 1._wp 371 en (ji,jj,itop) = z_en ; en (ji,jj,itopp1) = z_en 372 ENDIF 373 END_2D 374 ENDIF 329 375 ! 330 376 CASE ( 1 ) ! Neumann boundary condition (set d(e)/dz) … … 350 396 END_2D 351 397 ! 398 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 399 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 400 IF( mikt(ji,jj) > 1 )THEN 401 itop = mikt(ji,jj) ! k top w-point 402 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 403 ! ! mask at the 404 ! ocean surface 405 ! points 406 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 407 ! 408 ! Bottom level Dirichlet condition: 409 ! Bottom level (ibot) & Just above it (ibotm1) 410 ! Dirichlet ! Neumann 411 zd_lw(ji,jj,itop) = 0._wp ! ! Remove zd_up from zdiag 412 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) 413 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 414 en (ji,jj,itop) = z_en 415 ENDIF 416 END_2D 417 ENDIF 352 418 ! 353 419 END SELECT … … 605 671 END_2D 606 672 ! 673 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 674 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 675 IF ( mikt(ji,jj) > 1 ) THEN 676 itop = mikt(ji,jj) ! k top w-point 677 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 678 ! 679 zdep(ji,jj) = vkarmn * r_z0_top 680 psi (ji,jj,itop) = rc0**rpp * en(ji,jj,itop)**rmm *zdep(ji,jj)**rnn 681 zd_lw(ji,jj,itop) = 0._wp 682 zd_up(ji,jj,itop) = 0._wp 683 zdiag(ji,jj,itop) = 1._wp 684 ! 685 ! Just above last level, Dirichlet condition again (GOTM like) 686 zdep(ji,jj) = vkarmn * ( r_z0_top + e3t(ji,jj,itopp1,Kmm) ) 687 psi (ji,jj,itopp1) = rc0**rpp * en(ji,jj,itop )**rmm *zdep(ji,jj)**rnn 688 zd_lw(ji,jj,itopp1) = 0._wp 689 zd_up(ji,jj,itopp1) = 0._wp 690 zdiag(ji,jj,itopp1) = 1._wp 691 END IF 692 END_2D 693 END IF 694 ! 607 695 CASE ( 1 ) ! Neumman boundary condition 608 696 ! … … 629 717 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w(ji,jj,ibotm1,Kmm) 630 718 END_2D 719 ! 720 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 721 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 722 IF ( mikt(ji,jj) > 1 ) THEN 723 itop = mikt(ji,jj) ! k top w-point 724 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 725 ! 726 ! Bottom level Dirichlet condition: 727 zdep(ji,jj) = vkarmn * r_z0_top 728 psi (ji,jj,itop) = rc0**rpp * en(ji,jj,itop)**rmm *zdep(ji,jj)**rnn 729 ! 730 zd_lw(ji,jj,itop) = 0._wp 731 zd_up(ji,jj,itop) = 0._wp 732 zdiag(ji,jj,itop) = 1._wp 733 ! 734 ! Just below cavity level: Neumann condition with flux 735 ! injection 736 zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) ! Remove zd_up from zdiag 737 zd_up(ji,jj,itopp1) = 0._wp 738 ! 739 ! Set psi vertical flux below cavity: 740 zdep(ji,jj) = r_z0_top + 0.5_wp*e3t(ji,jj,itopp1,Kmm) 741 zflxb = rsbc_psi2 * ( p_avm(ji,jj,itop) + p_avm(ji,jj,itopp1)) & 742 & * (0.5_wp*(en(ji,jj,itop)+en(ji,jj,itopp1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 743 psi(ji,jj,itopp1) = psi(ji,jj,itopp1) + zflxb / e3w(ji,jj,itopp1,Kmm) 744 END IF 745 END_2D 746 END IF 747 631 748 ! 632 749 END SELECT … … 797 914 NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 798 915 & rn_clim_galp, ln_sigpsi, rn_hsro, rn_hsri, & 799 & rn_crban, rn_charn, rn_frac_hs,&916 & nn_mxlice, rn_crban, rn_charn, rn_frac_hs, & 800 917 & nn_bc_surf, nn_bc_bot, nn_z0_met, nn_z0_ice, & 801 918 & nn_stab_func, nn_clos … … 837 954 WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos 838 955 WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro 839 WRITE(numout,*) ' Ice-ocean roughness (used if nn_z0_ice/=0) rn_hsri = ', rn_hsri 956 WRITE(numout,*) ' type of scaling under sea-ice nn_mxlice = ', nn_mxlice 957 IF( nn_mxlice == 1 ) & 958 WRITE(numout,*) ' Ice-ocean roughness (used if nn_z0_ice/=0) rn_hsri = ', rn_hsri 959 SELECT CASE( nn_mxlice ) ! Type of scaling under sea-ice 960 CASE( 0 ) ; WRITE(numout,*) ' ==>>> No scaling under sea-ice' 961 CASE( 1 ) ; WRITE(numout,*) ' ==>>> scaling with constant sea-ice thickness' 962 CASE( 2 ) ; WRITE(numout,*) ' ==>>> scaling with mean sea-ice thickness' 963 CASE( 3 ) ; WRITE(numout,*) ' ==>>> scaling with max sea-ice thickness' 964 CASE DEFAULT 965 CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 ') 966 END SELECT 840 967 WRITE(numout,*) 841 968 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.