Changeset 7991 for branches/2017/dev_r7881_HPC09_ZDF
- Timestamp:
- 2017-05-01T16:21:42+02:00 (7 years ago)
- Location:
- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 1 added
- 4 edited
-
DIA/dia25h.F90 (modified) (6 diffs)
-
ZDF/zdfgls.F90 (modified) (26 diffs)
-
ZDF/zdfric.F90 (modified) (3 diffs)
-
ZDF/zdfsh2.F90 (added)
-
ZDF/zdftke.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90
r7990 r7991 9 9 USE dom_oce ! ocean space and time domain 10 10 USE zdf_oce ! ocean vertical physics 11 USE zdfgls , ONLY : hmx n11 USE zdfgls , ONLY : hmxl_n 12 12 USE in_out_manager ! I/O units 13 13 USE iom ! I/0 library … … 23 23 24 24 ! variables for calculating 25-hourly means 25 INTEGER , SAVE :: cnt_25h ! Counter for 25 hour means 25 INTEGER , SAVE :: cnt_25h ! Counter for 25 hour means 26 REAL(wp), SAVE :: r1_25 = 0.04_wp ! =1/25 26 27 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: tn_25h , sn_25h 27 28 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: sshn_25h … … 94 95 ! ------------------------- ! 95 96 cnt_25h = 1 ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 96 tn_25h (:,:,:) = tsb(:,:,:,jp_tem)97 sn_25h (:,:,:) = tsb(:,:,:,jp_sal)98 sshn_25h(:,:) = sshb(:,:)99 un_25h (:,:,:) = ub(:,:,:)100 vn_25h (:,:,:) = vb(:,:,:)101 wn_25h (:,:,:) = wn(:,:,:)102 avt_25h (:,:,:) = avt(:,:,:)103 avm_25h (:,:,:) = avm(:,:,:)97 tn_25h (:,:,:) = tsb (:,:,:,jp_tem) 98 sn_25h (:,:,:) = tsb (:,:,:,jp_sal) 99 sshn_25h(:,:) = sshb(:,:) 100 un_25h (:,:,:) = ub (:,:,:) 101 vn_25h (:,:,:) = vb (:,:,:) 102 wn_25h (:,:,:) = wn (:,:,:) 103 avt_25h (:,:,:) = avt (:,:,:) 104 avm_25h (:,:,:) = avm (:,:,:) 104 105 IF( ln_zdftke ) THEN 105 106 en_25h(:,:,:) = en(:,:,:) 106 107 ENDIF 107 108 IF( ln_zdfgls ) THEN 108 en_25h (:,:,:) = en(:,:,:)109 rmxln_25h(:,:,:) = hmx n(:,:,:)109 en_25h (:,:,:) = en (:,:,:) 110 rmxln_25h(:,:,:) = hmxl_n(:,:,:) 110 111 ENDIF 111 112 #if defined key_lim3 || defined key_lim2 … … 156 157 ENDIF 157 158 158 tn_25h (:,:,:) = tn_25h(:,:,:) + tsn(:,:,:,jp_tem)159 sn_25h (:,:,:) = sn_25h(:,:,:) + tsn(:,:,:,jp_sal)160 sshn_25h(:,:) = sshn_25h(:,:) + sshn(:,:)161 un_25h (:,:,:) = un_25h(:,:,:) + un(:,:,:)162 vn_25h (:,:,:) = vn_25h(:,:,:) + vn(:,:,:)163 wn_25h (:,:,:) = wn_25h(:,:,:) + wn(:,:,:)164 avt_25h (:,:,:) = avt_25h(:,:,:) + avt(:,:,:)165 avm_25h (:,:,:) = avm_25h(:,:,:) + avm(:,:,:)166 IF( ln_zdftke ) THEN 167 en_25h(:,:,:) = en_25h(:,:,:) + en(:,:,:)168 ENDIF 169 IF( ln_zdfgls ) THEN 170 en_25h (:,:,:) = en_25h(:,:,:) + en(:,:,:)171 rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + hmxn(:,:,:)159 tn_25h (:,:,:) = tn_25h (:,:,:) + tsn (:,:,:,jp_tem) 160 sn_25h (:,:,:) = sn_25h (:,:,:) + tsn (:,:,:,jp_sal) 161 sshn_25h(:,:) = sshn_25h(:,:) + sshn(:,:) 162 un_25h (:,:,:) = un_25h (:,:,:) + un (:,:,:) 163 vn_25h (:,:,:) = vn_25h (:,:,:) + vn (:,:,:) 164 wn_25h (:,:,:) = wn_25h (:,:,:) + wn (:,:,:) 165 avt_25h (:,:,:) = avt_25h (:,:,:) + avt (:,:,:) 166 avm_25h (:,:,:) = avm_25h (:,:,:) + avm (:,:,:) 167 IF( ln_zdftke ) THEN 168 en_25h(:,:,:) = en_25h (:,:,:) + en(:,:,:) 169 ENDIF 170 IF( ln_zdfgls ) THEN 171 en_25h (:,:,:) = en_25h (:,:,:) + en (:,:,:) 172 rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + hmxl_n(:,:,:) 172 173 ENDIF 173 174 cnt_25h = cnt_25h + 1 … … 187 188 ENDIF 188 189 ! 189 tn_25h (:,:,:) = tn_25h(:,:,:) / 25.0_wp190 sn_25h (:,:,:) = sn_25h(:,:,:) / 25.0_wp191 sshn_25h(:,:) = sshn_25h(:,:) / 25.0_wp192 un_25h (:,:,:) = un_25h(:,:,:) / 25.0_wp193 vn_25h (:,:,:) = vn_25h(:,:,:) / 25.0_wp194 wn_25h (:,:,:) = wn_25h(:,:,:) / 25.0_wp195 avt_25h (:,:,:) = avt_25h(:,:,:) / 25.0_wp196 avm_25h (:,:,:) = avm_25h(:,:,:) / 25.0_wp197 IF( ln_zdftke ) THEN 198 en_25h(:,:,:) = en_25h(:,:,:) / 25.0_wp199 ENDIF 200 IF( ln_zdfgls ) THEN 201 en_25h (:,:,:) = en_25h(:,:,:) / 25.0_wp202 rmxln_25h(:,:,:) = rmxln_25h(:,:,:) / 25.0_wp190 tn_25h (:,:,:) = tn_25h (:,:,:) * r1_25 191 sn_25h (:,:,:) = sn_25h (:,:,:) * r1_25 192 sshn_25h(:,:) = sshn_25h(:,:) * r1_25 193 un_25h (:,:,:) = un_25h (:,:,:) * r1_25 194 vn_25h (:,:,:) = vn_25h (:,:,:) * r1_25 195 wn_25h (:,:,:) = wn_25h (:,:,:) * r1_25 196 avt_25h (:,:,:) = avt_25h (:,:,:) * r1_25 197 avm_25h (:,:,:) = avm_25h (:,:,:) * r1_25 198 IF( ln_zdftke ) THEN 199 en_25h(:,:,:) = en_25h(:,:,:) * r1_25 200 ENDIF 201 IF( ln_zdfgls ) THEN 202 en_25h (:,:,:) = en_25h (:,:,:) * r1_25 203 rmxln_25h(:,:,:) = rmxln_25h(:,:,:) * r1_25 203 204 ENDIF 204 205 ! … … 236 237 ! 237 238 ! After the write reset the values to cnt=1 and sum values equal current value 238 tn_25h (:,:,:) = tsn(:,:,:,jp_tem)239 sn_25h (:,:,:) = tsn(:,:,:,jp_sal)240 sshn_25h(:,:) = sshn(:,:)241 un_25h (:,:,:) = un(:,:,:)242 vn_25h (:,:,:) = vn(:,:,:)243 wn_25h (:,:,:) = wn(:,:,:)244 avt_25h (:,:,:) = avt(:,:,:)245 avm_25h (:,:,:) = avm(:,:,:)239 tn_25h (:,:,:) = tsn (:,:,:,jp_tem) 240 sn_25h (:,:,:) = tsn (:,:,:,jp_sal) 241 sshn_25h(:,:) = sshn(:,:) 242 un_25h (:,:,:) = un (:,:,:) 243 vn_25h (:,:,:) = vn (:,:,:) 244 wn_25h (:,:,:) = wn (:,:,:) 245 avt_25h (:,:,:) = avt (:,:,:) 246 avm_25h (:,:,:) = avm (:,:,:) 246 247 IF( ln_zdftke ) THEN 247 248 en_25h(:,:,:) = en(:,:,:) 248 249 ENDIF 249 250 IF( ln_zdfgls ) THEN 250 en_25h (:,:,:) = en(:,:,:)251 rmxln_25h(:,:,:) = hmx n(:,:,:)251 en_25h (:,:,:) = en (:,:,:) 252 rmxln_25h(:,:,:) = hmxl_n(:,:,:) 252 253 ENDIF 253 254 cnt_25h = 1 -
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 -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
r7990 r7991 21 21 USE oce ! ocean dynamics and tracers variables 22 22 USE dom_oce ! ocean space and time domain variables 23 USE zdf_oce ! ocean vertical physics 23 USE zdf_oce ! vertical physics: variables 24 USE zdfsh2 ! vertical physics: shear production term of TKE 24 25 USE phycst ! physical constants 25 26 USE sbc_oce, ONLY : taum … … 104 105 INTEGER :: ji, jj, jk ! dummy loop indices 105 106 LOGICAL :: ll_av_weight = .TRUE. ! local logical 106 REAL(wp) :: zsh2, zcfRi, zav, zustar ! local scalars 107 REAL(wp), DIMENSION(jpi,jpj) :: zdu2, zdv2, zh_ekm ! 2D workspace 107 REAL(wp) :: zcfRi, zav, zustar ! local scalars 108 REAL(wp), DIMENSION(jpi,jpj) :: zh_ekm ! 2D workspace 109 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsh2 ! 3D - 108 110 !!---------------------------------------------------------------------- 109 111 ! 110 112 IF( nn_timing == 1 ) CALL timing_start('zdf_ric') 111 113 ! 112 DO jk = 2, jpkm1 !== avm and avt = F(Richardson number) ==! 113 ! 114 IF( ll_av_weight ) THEN !== viscosity weighted shear & stratification terms 115 ! 116 DO jj = 1, jpjm1 !* viscosity weighted shear term 114 ! !== avm and avt = F(Richardson number) ==! 115 ! 116 ! !* Shear production at uw- and vw-points (energy conserving form) 117 CALL zdf_sh2( zsh2 ) 118 ! 119 DO jk = 2, jpkm1 !* Vertical eddy viscosity and diffusivity coefficients 120 DO jj = 1, jpjm1 117 121 DO ji = 1, jpim1 118 zdu2(ji,jj) = 0.5 * ( avm(ji,jj,jk ) + avm(ji+1,jj,jk) ) * wumask(ji,jj,jk) & 119 & * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 120 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 121 zdv2(ji,jj) = 0.5 * ( avm(ji,jj,jk ) + avm(ji,jj+1,jk) ) * wumask(ji,jj,jk) & 122 & * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 123 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 124 END DO 125 END DO 126 DO jj = 2, jpjm1 127 DO ji = 2, jpim1 !* Richardson number dependent coefficient 128 ! ! shear of horizontal velocity 129 zsh2 = ( zdu2(ji-1,jj) + zdu2(ji,jj) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 130 & + ( zdv2(ji,jj-1) + zdv2(ji,jj) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 131 ! ! coefficient = F(richardson number) 132 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , avt(ji,jj,jk) * rn2(ji,jj,jk) / ( zsh2 + 1.e-20 ) ) ) 122 ! ! coefficient = F(richardson number) (avm-weighted Ri) 123 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( zsh2(ji,jj,jk) + 1.e-20 ) ) ) 124 zav = rn_avmri * zcfRi**nn_ric 133 125 ! 134 ! !* Vertical eddy viscosity and diffusivity coefficients 135 zav = rn_avmri * zcfRi**nn_ric 126 ! ! avm and avt coefficients 136 127 avm(ji,jj,jk) = MAX( zav , avmb(jk) ) * wmask(ji,jj,jk) 137 128 avt(ji,jj,jk) = MAX( zav * zcfRi , avtb(jk) ) * wmask(ji,jj,jk) 138 129 END DO 139 130 END DO 140 ELSE !== classical Richardson number definition ==!141 DO jj = 1, jpjm1 !* shear term142 DO ji = 1, jpim1143 zdu2(ji,jj) = 0.5 * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) &144 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )145 zdv2(ji,jj) = 0.5 * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) &146 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )147 END DO148 END DO149 DO jj = 2, jpjm1150 DO ji = 2, jpim1 !* Richardson number dependent coefficient151 ! ! shear of horizontal velocity152 zsh2 = ( zdu2(ji-1,jj) + zdu2(ji,jj) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &153 & + ( zdv2(ji,jj-1) + zdv2(ji,jj) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )154 ! ! coefficient = F(richardson number)155 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , rn2(ji,jj,jk) / ( zsh2 + 1.e-20 ) ) )156 !157 ! !* Vertical eddy viscosity and diffusivity coefficients158 zav = rn_avmri * zcfRi**nn_ric159 avm(ji,jj,jk) = MAX( zav , avmb(jk) ) * wmask(ji,jj,jk)160 avt(ji,jj,jk) = MAX( zav * zcfRi , avtb(jk) ) * wmask(ji,jj,jk)161 END DO162 END DO163 ENDIF164 !165 131 END DO 132 ! 133 !!gm BUG <<<<==== This param can't work at low latitude 134 !!gm it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! ) 166 135 ! 167 136 IF( ln_mldw ) THEN !== set a minimum value in the Ekman layer ==! … … 184 153 END DO 185 154 END DO 186 END 155 ENDIF 187 156 ! 188 157 IF( nn_timing == 1 ) CALL timing_stop('zdf_ric') -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r7990 r7991 43 43 USE sbc_oce ! surface boundary condition: ocean 44 44 USE zdf_oce ! vertical physics: ocean variables 45 USE zdfsh2 ! vertical physics: shear production term of TKE 45 46 USE zdfmxl ! vertical physics: mixed layer 46 47 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 89 90 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 90 91 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation 91 #if defined key_c1d92 ! !!** 1D cfg only ** ('key_c1d')93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales94 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers95 #endif96 92 97 93 !! * Substitutions … … 108 104 !! *** FUNCTION zdf_tke_alloc *** 109 105 !!---------------------------------------------------------------------- 110 ALLOCATE( & 111 #if defined key_c1d 112 & e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , & 113 & e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) , & 114 #endif 115 & htau (jpi,jpj) , dissl(jpi,jpj,jpk) , & 116 & apdlr(jpi,jpj,jpk) , STAT= zdf_tke_alloc ) 106 ALLOCATE( htau(jpi,jpj) , dissl(jpi,jpj,jpk) , apdlr(jpi,jpj,jpk) , STAT= zdf_tke_alloc ) 117 107 ! 118 108 IF( lk_mpp ) CALL mpp_sum ( zdf_tke_alloc ) … … 211 201 !! --------------------------------------------------------------------- 212 202 INTEGER :: ji, jj, jk ! dummy loop arguments 213 !!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! temporary scalar214 !!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! temporary scalar215 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 216 REAL(wp) :: z cdrag = 1.5e-3 ! drag coefficient217 REAL(wp) :: z bbrau, zesh2 ! temporary scalars218 REAL(wp) :: z fact1, zfact2, zfact3 ! - -219 REAL(wp) :: z tx2 , zty2 , zcof !- -220 REAL(wp) :: zt au , zdif !- -221 REAL(wp) :: z us , zwlc , zind !- -222 REAL(wp) :: z zd_up, zzd_lw !- -223 !!bfr REAL(wp) :: zebot !- -203 !!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! local integers 204 !!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! - - 205 !!bfr REAL(wp) :: zebot ! local scalars 206 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 207 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient 208 REAL(wp) :: zbbrau, zri ! local scalars 209 REAL(wp) :: zfact1, zfact2, zfact3 ! - - 210 REAL(wp) :: ztx2 , zty2 , zcof ! - - 211 REAL(wp) :: ztau , zdif ! - - 212 REAL(wp) :: zus , zwlc , zind ! - - 213 REAL(wp) :: zzd_up, zzd_lw ! - - 224 214 INTEGER , DIMENSION(jpi,jpj) :: imlc 225 215 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 226 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw , z3du, z3dv227 REAL(wp) :: zri ! local Richardson number216 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 217 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsh2 228 218 !!-------------------------------------------------------------------- 229 219 ! … … 231 221 ! 232 222 CALL wrk_alloc( jpi,jpj, zhlc ) 233 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw , z3du, z3dv)223 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 234 224 ! 235 225 zbbrau = rn_ebb / rau0 ! Local constant initialisation … … 328 318 ! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 329 319 ! 330 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 331 DO jj = 1, jpjm1 332 DO ji = 1, fs_jpim1 ! vector opt. 333 z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji+1,jj,jk) ) & 334 & * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 335 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) & 336 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 337 z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji,jj+1,jk) ) & 338 & * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 339 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) & 340 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 341 END DO 342 END DO 343 END DO 344 ! 345 IF( nn_pdl == 1 ) THEN !* Prandtl number case: compute apdlr 346 ! Note that zesh2 is also computed in the next loop. 347 ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops 320 ! !* Shear production at uw- and vw-points (energy conserving form) 321 CALL zdf_sh2( zsh2 ) 322 ! 323 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 348 324 DO jk = 2, jpkm1 349 325 DO jj = 2, jpjm1 350 DO ji = fs_2, fs_jpim1 ! vector opt. 351 ! ! shear prod. at w-point weightened by mask 352 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 353 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 354 ! ! local Richardson number 355 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear ) 326 DO ji = 2, jpim1 327 ! ! local Richardson number 328 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zsh2(ji,jj,jk) + rn_bshear ) 329 ! ! inverse of Prandtl number 356 330 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) 357 358 END DO 359 END DO 360 END DO 361 ! 331 END DO 332 END DO 333 END DO 362 334 ENDIF 363 335 ! … … 373 345 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 374 346 ! 375 ! ! shear prod. at w-point weightened by mask376 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &377 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )378 !379 347 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 380 348 zd_lw(ji,jj,jk) = zzd_lw 381 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)349 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 382 350 ! 383 351 ! ! right hand side in en 384 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 385 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 386 & * wmask(ji,jj,jk) 352 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zsh2(ji,jj,jk) & ! shear 353 & - avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 354 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 355 & ) * wmask(ji,jj,jk) 387 356 END DO 388 357 END DO … … 470 439 ! 471 440 CALL wrk_dealloc( jpi,jpj, zhlc ) 472 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw , z3du, z3dv)441 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 473 442 ! 474 443 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') … … 631 600 END SELECT 632 601 ! 633 # if defined key_c1d634 e_dis(:,:,:) = zmxld(:,:,:) ! c1d configuration : save mixing and dissipation turbulent length scales635 e_mix(:,:,:) = zmxlm(:,:,:)636 # endif637 602 638 603 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 657 622 DO ji = fs_2, fs_jpim1 ! vector opt. 658 623 avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 659 # if defined key_c1d660 e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number661 # endif662 624 END DO 663 625 END DO
Note: See TracChangeset
for help on using the changeset viewer.
