- Timestamp:
- 2017-05-01T16:21:42+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r7990 r7991 19 19 USE domvvl ! ocean space and time domain : variable volume layer 20 20 USE zdf_oce ! ocean vertical physics 21 USE zdfsh2 ! vertical physics: shear production term of TKE 21 22 USE zdfbfr ! bottom friction (only for rn_bfrz0) 22 23 USE sbc_oce ! surface boundary condition: ocean … … 42 43 43 44 ! 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hmx n !: now mixing length45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hmxl_n !: now mixing length 45 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zwall !: wall function 46 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points … … 114 115 !! *** FUNCTION zdf_gls_alloc *** 115 116 !!---------------------------------------------------------------------- 116 ALLOCATE( hmx n(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,&117 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) ,STAT= zdf_gls_alloc )117 ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustars2(jpi,jpj) , & 118 & zwall (jpi,jpj,jpk) , ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 118 119 ! 119 120 IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc ) … … 141 142 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) 142 143 REAL(wp), POINTER, DIMENSION(:,:,:) :: eb ! tke at time before 143 REAL(wp), POINTER, DIMENSION(:,:,:) :: mxlb! mixing length at time before144 REAL(wp), POINTER, DIMENSION(:,:,:) :: hmxl_b ! mixing length at time before 144 145 REAL(wp), POINTER, DIMENSION(:,:,:) :: shear ! vertical shear 145 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: eps ! dissipation rate … … 149 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_b ! element of the second matrix diagonal 150 151 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_c ! element of the third matrix diagonal 151 REAL(wp), POINTER, DIMENSION(:,:,:) :: z 3du, z3dv ! u- and v-shears152 REAL(wp), POINTER, DIMENSION(:,:,:) :: zstt, zstm ! stability function on tracer and momentum 152 153 !!-------------------------------------------------------------------- 153 154 ! … … 155 156 ! 156 157 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 157 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )158 CALL wrk_alloc( jpi,jpj,jpk, z 3du, z3dv)158 CALL wrk_alloc( jpi,jpj,jpk, eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 159 CALL wrk_alloc( jpi,jpj,jpk, zstt, zstm ) 159 160 160 161 ! Preliminary computing … … 198 199 END SELECT 199 200 200 DO jk = 2, jpkm1 !== Compute shear and dissipation rate ==! 201 ! !== Compute shear and dissipation rate ==! 202 CALL zdf_sh2( shear ) 203 204 205 DO jk = 2, jpkm1 !== Compute dissipation rate ==! 201 206 DO jj = 1, jpjm1 202 207 DO ji = 1, jpim1 203 z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji+1,jj,jk) ) & 204 & * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 205 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 206 z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji,jj+1,jk) ) & 207 & * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 208 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 209 ! 210 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / hmxn(ji,jj,jk) 208 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 211 209 END DO 212 210 END DO … … 214 212 215 213 ! Save tke at before time step 216 eb (:,:,:) = en(:,:,:)217 mxlb(:,:,:) = hmxn(:,:,:)214 eb (:,:,:) = en (:,:,:) 215 hmxl_b(:,:,:) = hmxl_n(:,:,:) 218 216 219 217 IF( nn_clos == 0 ) THEN ! Mellor-Yamada … … 221 219 DO jj = 2, jpjm1 222 220 DO ji = fs_2, fs_jpim1 ! vector opt. 223 zup = hmx n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)221 zup = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 224 222 zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 225 223 zcoef = ( zup / MAX( zdown, rsmall ) ) … … 247 245 DO ji = 2, jpim1 248 246 ! 249 ! shear prod. at w-point weightened by mask 250 shear(ji,jj,jk) = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 251 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 252 ! 253 ! stratif. destruction 254 buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk) 255 ! 256 ! shear prod. - stratif. destruction 257 diss = eps(ji,jj,jk) 247 buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction 248 ! 249 diss = eps(ji,jj,jk) ! dissipation 258 250 ! 259 251 dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy ) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) … … 344 336 ! 345 337 CASE ( 0 ) ! Dirichlet 346 ! ! en(ibot) = u*^2 / Co2 and hmx n(ibot) = rn_lmin338 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 347 339 ! ! Balance between the production and the dissipation terms 348 340 DO jj = 2, jpjm1 … … 425 417 DO jj = 2, jpjm1 426 418 DO ji = fs_2, fs_jpim1 ! vector opt. 427 psi(ji,jj,jk) = eb(ji,jj,jk) * mxlb(ji,jj,jk)419 psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 428 420 END DO 429 421 END DO … … 443 435 DO jj = 2, jpjm1 444 436 DO ji = fs_2, fs_jpim1 ! vector opt. 445 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * mxlb(ji,jj,jk) )437 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 446 438 END DO 447 439 END DO … … 452 444 DO jj = 2, jpjm1 453 445 DO ji = fs_2, fs_jpim1 ! vector opt. 454 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn446 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 455 447 END DO 456 448 END DO … … 563 555 ! 564 556 CASE ( 0 ) ! Dirichlet 565 ! ! en(ibot) = u*^2 / Co2 and hmx n(ibot) = vkarmn * rn_bfrz0557 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * rn_bfrz0 566 558 ! ! Balance between the production and the dissipation terms 567 559 DO jj = 2, jpjm1 … … 686 678 ! Limit dissipation rate under stable stratification 687 679 ! -------------------------------------------------- 688 DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmx n at the same time680 DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time 689 681 DO jj = 2, jpjm1 690 682 DO ji = fs_2, fs_jpim1 ! vector opt. 691 683 ! limitation 692 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin )693 hmx n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)684 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 685 hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 694 686 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 695 687 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 696 IF (ln_length_lim) hmxn(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxn(ji,jj,jk) )688 IF( ln_length_lim ) hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 697 689 END DO 698 690 END DO … … 710 702 DO ji = fs_2, fs_jpim1 ! vector opt. 711 703 ! zcof = l²/q² 712 zcof = mxlb(ji,jj,jk) * mxlb(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )704 zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 713 705 ! Gh = -N²l²/q² 714 706 gh = - rn2(ji,jj,jk) * zcof … … 719 711 sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh) 720 712 ! 721 ! Store stability function in z 3du and z3dv722 z 3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)723 z 3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)713 ! Store stability function in zstt and zstm 714 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 715 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 724 716 END DO 725 717 END DO … … 731 723 DO ji = fs_2, fs_jpim1 ! vector opt. 732 724 ! zcof = l²/q² 733 zcof = mxlb(ji,jj,jk)*mxlb(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )725 zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 734 726 ! Gh = -N²l²/q² 735 727 gh = - rn2(ji,jj,jk) * zcof … … 747 739 sh = (rs4 - rs5*gh + rs6*gm) / rcff 748 740 ! 749 ! Store stability function in z 3du and z3dv750 z 3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)751 z 3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)741 ! Store stability function in zstt and zstm 742 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 743 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 752 744 END DO 753 745 END DO … … 759 751 ! Lines below are useless if GOTM style Dirichlet conditions are used 760 752 761 z3dv(:,:,1) = z3dv(:,:,2) 762 753 zstm(:,:,1) = zstm(:,:,2) 754 755 !!gm should be done for ISF (top boundary cond.) 763 756 DO jj = 2, jpjm1 764 757 DO ji = fs_2, fs_jpim1 ! vector opt. 765 z 3dv(ji,jj,mbkt(ji,jj)+1) = z3dv(ji,jj,mbkt(ji,jj))758 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 766 759 END DO 767 760 END DO … … 772 765 DO jj = 2, jpjm1 773 766 DO ji = fs_2, fs_jpim1 ! vector opt. 774 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmx n(ji,jj,jk)775 zav = zsqen * z 3du(ji,jj,jk)767 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 768 zav = zsqen * zstt(ji,jj,jk) 776 769 avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 777 zav = zsqen * z 3dv(ji,jj,jk)770 zav = zsqen * zstm(ji,jj,jk) 778 771 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 779 772 END DO … … 795 788 ! 796 789 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 797 CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )798 CALL wrk_dealloc( jpi,jpj,jpk, z 3du, z3dv)790 CALL wrk_dealloc( jpi,jpj,jpk, eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 791 CALL wrk_dealloc( jpi,jpj,jpk, zstt, zstm ) 799 792 ! 800 793 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls') … … 1100 1093 1101 1094 ! !* read or initialize all required files 1102 CALL gls_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, hmx n)1095 CALL gls_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, hmxl_n) 1103 1096 ! 1104 1097 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls_init') … … 1129 1122 ! ! --------------- 1130 1123 IF( ln_rstart ) THEN !* Read the restart file 1131 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. )1132 id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. )1133 id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. )1134 id4 = iom_varid( numror, 'hmx n', ldstop = .FALSE. )1124 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 1125 id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. ) 1126 id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. ) 1127 id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. ) 1135 1128 ! 1136 1129 IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! all required arrays exist 1137 CALL iom_get( numror, jpdom_autoglo, 'en' , en)1138 CALL iom_get( numror, jpdom_autoglo, 'avt_k' , avt_k)1139 CALL iom_get( numror, jpdom_autoglo, 'avm_k' , avm_k)1140 CALL iom_get( numror, jpdom_autoglo, 'hmx n' , hmxn)1130 CALL iom_get( numror, jpdom_autoglo, 'en' , en ) 1131 CALL iom_get( numror, jpdom_autoglo, 'avt_k' , avt_k ) 1132 CALL iom_get( numror, jpdom_autoglo, 'avm_k' , avm_k ) 1133 CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n ) 1141 1134 ELSE 1142 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmx n computed by iterative loop'1143 en (:,:,:) = rn_emin1144 hmx n(:,:,:) = 0.05_wp1145 avt_k (:,:,:) = avt(:,:,:)1146 avm_k (:,:,:) = avm(:,:,:)1135 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxl_n computed by iterative loop' 1136 en (:,:,:) = rn_emin 1137 hmxl_n(:,:,:) = 0.05_wp 1138 avt_k (:,:,:) = avt(:,:,:) 1139 avm_k (:,:,:) = avm(:,:,:) 1147 1140 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_gls( jit ) ; END DO 1148 1141 ENDIF 1149 1142 ELSE !* Start from rest 1150 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmx n by background values'1151 en (:,:,:) = rn_emin1152 hmx n(:,:,:) = 0.05_wp1143 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxl_n by background values' 1144 en (:,:,:) = rn_emin 1145 hmxl_n(:,:,:) = 0.05_wp 1153 1146 DO jk = 1, jpk 1154 1147 avt_k(:,:,jk) = avtb(jk) * wmask(:,:,jk) … … 1160 1153 ! ! ------------------- 1161 1154 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1162 CALL iom_rstput( kt, nitrst, numrow, 'en' , en)1163 CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k)1164 CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k)1165 CALL iom_rstput( kt, nitrst, numrow, 'hmx n' , hmxn)1155 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1156 CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k ) 1157 CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k ) 1158 CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n ) 1166 1159 ! 1167 1160 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
