- Timestamp:
- 2017-04-29T17:24:54+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
r7953 r7990 5 5 !! turbulent closure parameterization 6 6 !!====================================================================== 7 !! History : 3.0 ! 2009-09 (G. Reffray) Original code 8 !! 3.3 ! 2010-10 (C. Bricaud) Add in the reference 7 !! History : 3.0 ! 2009-09 (G. Reffray) Original code 8 !! 3.3 ! 2010-10 (C. Bricaud) Add in the reference 9 !! 4.0 ! 2017-04 (G. Madec) remove CPP keys & avm at t-point only 9 10 !!---------------------------------------------------------------------- 10 11 … … 36 37 PRIVATE 37 38 38 PUBLIC zdf_gls ! routine called in step module39 PUBLIC zdf_gls_init ! routine called in zdfphy module40 PUBLIC gls_rst ! routine called in zdfphy module39 PUBLIC zdf_gls ! called in zdfphy 40 PUBLIC zdf_gls_init ! called in zdfphy 41 PUBLIC gls_rst ! called in zdfphy 41 42 42 43 ! 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mxln !: now mixing length44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hmxn !: now mixing length 44 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zwall !: wall function 45 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points … … 113 114 !! *** FUNCTION zdf_gls_alloc *** 114 115 !!---------------------------------------------------------------------- 115 ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , &116 ALLOCATE( hmxn(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , & 116 117 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 117 118 ! … … 148 149 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_b ! element of the second matrix diagonal 149 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_c ! element of the third matrix diagonal 151 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3du, z3dv ! u- and v-shears 150 152 !!-------------------------------------------------------------------- 151 153 ! 152 IF( nn_timing == 1 ) CALL timing_start('zdf_gls')154 IF( nn_timing == 1 ) CALL timing_start('zdf_gls') 153 155 ! 154 156 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 155 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 156 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, z3du, z3dv ) 159 157 160 ! Preliminary computing 158 161 159 162 ustars2 = 0._wp ; ustarb2 = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp 160 163 161 IF( kt /= nit000 ) THEN ! restore before value to compute tke 162 avt (:,:,:) = avt_k (:,:,:) 163 avm (:,:,:) = avm_k (:,:,:) 164 avmu(:,:,:) = avmu_k(:,:,:) 165 avmv(:,:,:) = avmv_k(:,:,:) 166 ENDIF 164 ! restore before values computed GLS alone 165 avt(:,:,:) = avt_k(:,:,:) 166 avm(:,:,:) = avm_k(:,:,:) 167 167 168 168 ! Compute surface and bottom friction at T-points … … 183 183 END DO 184 184 185 ! Set surface roughness length 186 SELECT CASE ( nn_z0_met ) 187 ! 188 CASE ( 0 ) ! Constant roughness 185 SELECT CASE ( nn_z0_met ) !== Set surface roughness length ==! 186 CASE ( 0 ) ! Constant roughness 189 187 zhsro(:,:) = rn_hsro 190 188 CASE ( 1 ) ! Standard Charnock formula 191 189 zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro) 192 190 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 191 !!gm zcof = 2._wp * 0.6_wp / 28._wp 192 !!gm zdep(:,:) = 30._wp * TANH( zcof/ SQRT( MAX(ustars2(:,:),rsmall) ) ) ! Wave age (eq. 10) 193 193 zdep(:,:) = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall)))) ! Wave age (eq. 10) 194 194 zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 195 195 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 196 !!gm BUG missing a multiplicative coefficient.... 196 197 zhsro(:,:) = hsw(:,:) 197 198 END SELECT 198 199 199 ! Compute shear and dissipation rate200 DO jk = 2, jpkm1201 DO jj = 2, jpjm1202 DO ji = fs_2, fs_jpim1 ! vector opt.203 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) &204 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) &205 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )206 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) &207 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) &208 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )209 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)200 DO jk = 2, jpkm1 !== Compute shear and dissipation rate ==! 201 DO jj = 1, jpjm1 202 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) 210 211 END DO 211 212 END DO 212 213 END DO 213 !214 ! Lateral boundary conditions (avmu,avmv) (sign unchanged)215 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. )216 214 217 215 ! Save tke at before time step 218 216 eb (:,:,:) = en (:,:,:) 219 mxlb(:,:,:) = mxln(:,:,:)217 mxlb(:,:,:) = hmxn(:,:,:) 220 218 221 219 IF( nn_clos == 0 ) THEN ! Mellor-Yamada … … 223 221 DO jj = 2, jpjm1 224 222 DO ji = fs_2, fs_jpim1 ! vector opt. 225 zup = mxln(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)223 zup = hmxn(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 226 224 zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 227 225 zcoef = ( zup / MAX( zdown, rsmall ) ) … … 247 245 DO jk = 2, jpkm1 248 246 DO jj = 2, jpjm1 249 DO ji = fs_2, fs_jpim1 ! vector opt.247 DO ji = 2, jpim1 250 248 ! 251 249 ! shear prod. at w-point weightened by mask 252 shear(ji,jj,jk) = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &253 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )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) ) 254 252 ! 255 253 ! stratif. destruction … … 278 276 ! building the matrix 279 277 zcof = rfact_tke * tmask(ji,jj,jk) 280 ! 281 ! lower diagonal 282 z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 283 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 284 ! 285 ! upper diagonal 286 z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 287 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 288 ! 289 ! diagonal 290 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & 291 & + rdt * zdiss * tmask(ji,jj,jk) 292 ! 293 ! right hand side in en 294 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk) 278 ! ! lower diagonal 279 z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk ) + avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 280 ! ! upper diagonal 281 z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) + avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 282 ! ! diagonal 283 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 284 ! ! right hand side in en 285 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 295 286 END DO 296 287 END DO … … 300 291 ! 301 292 ! Set surface condition on zwall_psi (1 at the bottom) 302 zwall_psi(:,:, 1) = zwall_psi(:,:,2)303 zwall_psi(:,:,jpk) = 1. 293 zwall_psi(:,:, 1 ) = zwall_psi(:,:,2) 294 zwall_psi(:,:,jpk) = 1._wp 304 295 ! 305 296 ! Surface boundary condition on tke … … 353 344 ! 354 345 CASE ( 0 ) ! Dirichlet 355 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin346 ! ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = rn_lmin 356 347 ! ! Balance between the production and the dissipation terms 357 348 DO jj = 2, jpjm1 … … 503 494 ! building the matrix 504 495 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 505 ! lower diagonal 506 z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 507 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 508 ! upper diagonal 509 z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 510 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 511 ! diagonal 512 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & 513 & + rdt * zdiss * tmask(ji,jj,jk) 514 ! 515 ! right hand side in psi 516 psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk) 496 ! ! lower diagonal 497 z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk ) + avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 498 ! ! upper diagonal 499 z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) + avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 500 ! ! diagonal 501 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 502 ! ! right hand side in psi 503 psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 517 504 END DO 518 505 END DO … … 565 552 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 566 553 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 567 568 554 ! 569 555 ! … … 577 563 ! 578 564 CASE ( 0 ) ! Dirichlet 579 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0565 ! ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = vkarmn * rn_bfrz0 580 566 ! ! Balance between the production and the dissipation terms 581 567 DO jj = 2, jpjm1 … … 700 686 ! Limit dissipation rate under stable stratification 701 687 ! -------------------------------------------------- 702 DO jk = 1, jpkm1 ! Note that this set boundary conditions on mxln at the same time688 DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxn at the same time 703 689 DO jj = 2, jpjm1 704 690 DO ji = fs_2, fs_jpim1 ! vector opt. 705 691 ! limitation 706 692 eps(ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 707 mxln(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)693 hmxn(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 708 694 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 709 695 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 710 IF (ln_length_lim) mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) )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) ) 711 697 END DO 712 698 END DO … … 733 719 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) 734 720 ! 735 ! Store stability function in avmu and avmv736 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)737 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)721 ! Store stability function in z3du and z3dv 722 z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 723 z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 738 724 END DO 739 725 END DO … … 761 747 sh = (rs4 - rs5*gh + rs6*gm) / rcff 762 748 ! 763 ! Store stability function in avmu and avmv764 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)765 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)749 ! Store stability function in z3du and z3dv 750 z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 751 z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 766 752 END DO 767 753 END DO … … 773 759 ! Lines below are useless if GOTM style Dirichlet conditions are used 774 760 775 avmv(:,:,1) = avmv(:,:,2)761 z3dv(:,:,1) = z3dv(:,:,2) 776 762 777 763 DO jj = 2, jpjm1 778 764 DO ji = fs_2, fs_jpim1 ! vector opt. 779 avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj))765 z3dv(ji,jj,mbkt(ji,jj)+1) = z3dv(ji,jj,mbkt(ji,jj)) 780 766 END DO 781 767 END DO … … 786 772 DO jj = 2, jpjm1 787 773 DO ji = fs_2, fs_jpim1 ! vector opt. 788 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * mxln(ji,jj,jk)789 zav = zsqen * avmu(ji,jj,jk)790 avt(ji,jj,jk) = MAX( zav, avtb(jk) ) *tmask(ji,jj,jk) ! apply mask for zdfmxl routine791 zav = zsqen * avmv(ji,jj,jk)792 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom774 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxn(ji,jj,jk) 775 zav = zsqen * z3du(ji,jj,jk) 776 avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 777 zav = zsqen * z3dv(ji,jj,jk) 778 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 793 779 END DO 794 780 END DO 795 781 END DO 796 !797 ! Lateral boundary conditions (sign unchanged)798 782 avt(:,:,1) = 0._wp 783 !!gm I'm sure this is needed to compute the shear term at next time-step 799 784 CALL lbc_lnk( avm, 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. ) 800 801 DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points 802 DO jj = 2, jpjm1 803 DO ji = fs_2, fs_jpim1 ! vector opt. 804 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 805 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 806 END DO 807 END DO 808 END DO 809 avmu(:,:,1) = 0._wp ; avmv(:,:,1) = 0._wp ! set surface to zero 810 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 811 785 ! 812 786 IF(ln_ctl) THEN 813 CALL prt_ctl( tab3d_1=en , clinfo1=' gls - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 814 CALL prt_ctl( tab3d_1=avmu, clinfo1=' gls - u: ', mask1=umask, & 815 & tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk ) 787 CALL prt_ctl( tab3d_1=en , clinfo1=' gls - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 788 CALL prt_ctl( tab3d_1=avm, clinfo1=' gls - m: ', ovlap=1, kdim=jpk ) 816 789 ENDIF 817 790 ! 818 avt_k (:,:,:) = avt (:,:,:)819 avm_k (:,:,:) = avm(:,:,:)820 avmu_k(:,:,:) = avmu(:,:,:)821 avmv_k(:,:,:) = avmv(:,:,:)791 avt_k(:,:,:) = avt(:,:,:) !== store avt, avm values computed by GLS only ==! 792 avm_k(:,:,:) = avm(:,:,:) 793 ! 794 IF( lrst_oce ) CALL gls_rst( kt, 'WRITE' ) ! write the GLS info in the restart file 822 795 ! 823 796 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 824 797 CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 825 !826 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls')827 !798 CALL wrk_dealloc( jpi,jpj,jpk, z3du, z3dv ) 799 ! 800 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls') 828 801 ! 829 802 END SUBROUTINE zdf_gls … … 835 808 !! 836 809 !! ** Purpose : Initialization of the vertical eddy diffivity and 837 !! viscosity when using a glsturbulent closure scheme810 !! viscosity computed using a GLS turbulent closure scheme 838 811 !! 839 812 !! ** Method : Read the namzdf_gls namelist and check the parameters 840 !! called at the first timestep (nit000)841 813 !! 842 814 !! ** input : Namlist namzdf_gls … … 892 864 ENDIF 893 865 894 ! !* allocate glsarrays866 ! !* allocate GLS arrays 895 867 IF( zdf_gls_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' ) 896 868 … … 1120 1092 rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 1121 1093 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1122 1094 ! 1123 1095 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1124 1096 rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1125 1097 ! 1126 1098 ! !* Wall proximity function 1127 1099 zwall (:,:,:) = 1._wp * tmask(:,:,:) 1128 1100 1129 ! !* set vertical eddy coef. to the background value 1130 DO jk = 1, jpk 1131 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 1132 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 1133 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 1134 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 1135 END DO 1136 ! 1137 CALL gls_rst( nit000, 'READ' ) !* read or initialize all required files 1101 ! !* read or initialize all required files 1102 CALL gls_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, hmxn) 1138 1103 ! 1139 1104 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls_init') … … 1152 1117 !! set to rn_emin or recomputed (nn_igls/=0) 1153 1118 !!---------------------------------------------------------------------- 1154 INTEGER , INTENT(in) :: kt 1155 CHARACTER(len=*), INTENT(in) :: cdrw 1119 INTEGER , INTENT(in) :: kt ! ocean time-step 1120 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1156 1121 ! 1157 1122 INTEGER :: jit, jk ! dummy loop indices 1158 INTEGER :: id1, id2, id3, id4 , id5, id61123 INTEGER :: id1, id2, id3, id4 1159 1124 INTEGER :: ji, jj, ikbu, ikbv 1160 1125 REAL(wp):: cbx, cby … … 1165 1130 IF( ln_rstart ) THEN !* Read the restart file 1166 1131 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 1167 id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. ) 1168 id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. ) 1169 id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. ) 1170 id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. ) 1171 id6 = iom_varid( numror, 'mxln' , 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, 'hmxn' , ldstop = .FALSE. ) 1172 1135 ! 1173 IF( MIN( id1, id2, id3, id4, id5, id6 ) > 0 ) THEN ! all required arrays exist 1174 CALL iom_get( numror, jpdom_autoglo, 'en' , en ) 1175 CALL iom_get( numror, jpdom_autoglo, 'avt' , avt ) 1176 CALL iom_get( numror, jpdom_autoglo, 'avm' , avm ) 1177 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu ) 1178 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv ) 1179 CALL iom_get( numror, jpdom_autoglo, 'mxln' , mxln ) 1136 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, 'hmxn' , hmxn ) 1180 1141 ELSE 1181 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 1182 en (:,:,:) = rn_emin 1183 mxln(:,:,:) = 0.05 1184 avt_k (:,:,:) = avt (:,:,:) 1185 avm_k (:,:,:) = avm (:,:,:) 1186 avmu_k(:,:,:) = avmu(:,:,:) 1187 avmv_k(:,:,:) = avmv(:,:,:) 1142 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxn computed by iterative loop' 1143 en (:,:,:) = rn_emin 1144 hmxn (:,:,:) = 0.05_wp 1145 avt_k(:,:,:) = avt(:,:,:) 1146 avm_k(:,:,:) = avm(:,:,:) 1188 1147 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_gls( jit ) ; END DO 1189 1148 ENDIF 1190 1149 ELSE !* Start from rest 1191 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'1150 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxn by background values' 1192 1151 en (:,:,:) = rn_emin 1193 mxln(:,:,:) = 0.05 1152 hmxn(:,:,:) = 0.05_wp 1153 DO jk = 1, jpk 1154 avt_k(:,:,jk) = avtb(jk) * wmask(:,:,jk) 1155 avm_k(:,:,jk) = avmb(jk) * wmask(:,:,jk) 1156 END DO 1194 1157 ENDIF 1195 1158 ! … … 1197 1160 ! ! ------------------- 1198 1161 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1199 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1200 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 1201 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 1202 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1203 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 1204 CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln ) 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, 'hmxn' , hmxn ) 1205 1166 ! 1206 1167 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
