- Timestamp:
- 2017-06-25T12:26:32+02:00 (7 years ago)
- Location:
- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 2 deleted
- 6 edited
-
dynbfr.F90 (modified) (5 diffs)
-
dynhpg.F90 (modified) (1 diff)
-
dynldf.F90 (modified) (1 diff)
-
dynldf_iso.F90 (modified) (6 diffs)
-
dynspg_ts.F90 (modified) (29 diffs)
-
dynzdf.F90 (modified) (5 diffs)
-
dynzdf_exp.F90 (deleted)
-
dynzdf_imp.F90 (deleted)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
r7753 r8215 5 5 !!============================================================================== 6 6 !! History : 3.2 ! 2008-11 (A. C. Coward) Original code 7 !! 3.4 ! 2011-09 (H. Liu) Make it consistent with semi-implicit 8 !! Bottom friction (ln_bfrimp = .true.)7 !! 3.4 ! 2011-09 (H. Liu) Make it consistent with semi-implicit Bottom friction (ln_drgimp =T) 8 !! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90) 9 9 !!---------------------------------------------------------------------- 10 10 … … 14 14 USE oce ! ocean dynamics and tracers variables 15 15 USE dom_oce ! ocean space and time domain variables 16 USE zdf_oce ! ocean vertical physicsvariables17 USE zdf bfr ! ocean bottom friction variables16 USE zdf_oce ! vertical physics: variables 17 USE zdfdrg ! vertical physics: top/bottom drag coef. 18 18 USE trd_oce ! trends: ocean variables 19 19 USE trddyn ! trend manager: dynamics 20 ! 20 21 USE in_out_manager ! I/O manager 21 22 USE prtctl ! Print control 22 23 USE timing ! Timing 23 USE wrk_nemo ! Memory Allocation24 24 25 25 IMPLICIT NONE … … 31 31 # include "vectopt_loop_substitute.h90" 32 32 !!---------------------------------------------------------------------- 33 !! NEMO/OPA 3.3 , NEMO Consortium (2010)33 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 34 34 !! $Id$ 35 35 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 43 43 !! ** Purpose : compute the bottom friction ocean dynamics physics. 44 44 !! 45 !! only for explicit bottom friction form 46 !! implicit bfr is implemented in dynzdf_imp 47 !! 45 48 !! ** Action : (ua,va) momentum trend increased by bottom friction trend 46 49 !!--------------------------------------------------------------------- … … 50 53 INTEGER :: ikbu, ikbv ! local integers 51 54 REAL(wp) :: zm1_2dt ! local scalar 52 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 55 REAL(wp) :: zCdu, zCdv ! - - 56 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv 53 57 !!--------------------------------------------------------------------- 54 58 ! 55 59 IF( nn_timing == 1 ) CALL timing_start('dyn_bfr') 56 60 ! 57 !!gm issue: better to put the logical in step to control the call of zdf_bfr 58 !! ==> change the logical from ln_bfrimp to ln_bfr_exp !! 59 IF( .NOT.ln_bfrimp) THEN ! only for explicit bottom friction form 60 ! implicit bfr is implemented in dynzdf_imp 61 !!gm bug : time step is only rdt (not 2 rdt if euler start !) 62 zm1_2dt = - 1._wp / ( 2._wp * rdt ) 61 63 62 !!gm bug : time step is only rdt (not 2 rdt if euler start !) 63 zm1_2dt = - 1._wp / ( 2._wp * rdt ) 64 65 IF( l_trddyn ) THEN ! trends: store the input trends 66 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 67 ztrdu(:,:,:) = ua(:,:,:) 68 ztrdv(:,:,:) = va(:,:,:) 69 ENDIF 64 IF( l_trddyn ) THEN ! trends: store the input trends 65 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 66 ztrdu(:,:,:) = ua(:,:,:) 67 ztrdv(:,:,:) = va(:,:,:) 68 ENDIF 70 69 71 70 72 DO jj = 2, jpjm1 73 DO ji = 2, jpim1 74 ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels 75 ikbv = mbkv(ji,jj) 76 ! 77 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 78 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( bfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) 79 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( bfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) 71 DO jj = 2, jpjm1 72 DO ji = 2, jpim1 73 ikbu = mbku(ji,jj) ! deepest wet ocean u- & v-levels 74 ikbv = mbkv(ji,jj) 75 ! 76 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 77 zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u_n(ji,jj,ikbu) 78 zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_n(ji,jj,ikbv) 79 ! 80 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( zCdu , zm1_2dt ) * ub(ji,jj,ikbu) 81 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( zCdv , zm1_2dt ) * vb(ji,jj,ikbv) 82 END DO 83 END DO 84 ! 85 IF( ln_isfcav ) THEN ! ocean cavities 86 DO jj = 2, jpjm1 87 DO ji = 2, jpim1 88 ikbu = miku(ji,jj) ! first wet ocean u- & v-levels 89 ikbv = mikv(ji,jj) 90 ! 91 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 92 zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u_n(ji,jj,ikbu) ! NB: Cdtop masked 93 zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v_n(ji,jj,ikbv) 94 ! 95 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( zCdu , zm1_2dt ) * ub(ji,jj,ikbu) 96 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( zCdv , zm1_2dt ) * vb(ji,jj,ikbv) 80 97 END DO 81 END DO 82 ! 83 IF( ln_isfcav ) THEN ! ocean cavities 84 DO jj = 2, jpjm1 85 DO ji = 2, jpim1 86 ! (ISF) stability criteria for top friction 87 ikbu = miku(ji,jj) ! first wet ocean u- & v-levels 88 ikbv = mikv(ji,jj) 89 ! 90 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 91 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( tfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) & 92 & * (1.-umask(ji,jj,1)) 93 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( tfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) & 94 & * (1.-vmask(ji,jj,1)) 95 ! (ISF) 96 END DO 97 END DO 98 END IF 99 ! 100 IF( l_trddyn ) THEN ! trends: send trends to trddyn for further diagnostics 101 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 102 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 103 CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt ) 104 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 105 ENDIF 106 ! ! print mean trends (used for debugging) 107 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr - Ua: ', mask1=umask, & 108 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 109 ! 110 ENDIF ! end explicit bottom friction 98 END DO 99 ENDIF 100 ! 101 IF( l_trddyn ) THEN ! trends: send trends to trddyn for further diagnostics 102 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 103 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 104 CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt ) 105 DEALLOCATE( ztrdu, ztrdv ) 106 ENDIF 107 ! ! print mean trends (used for debugging) 108 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr - Ua: ', mask1=umask, & 109 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 111 110 ! 112 111 IF( nn_timing == 1 ) CALL timing_stop('dyn_bfr') -
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r7761 r8215 1457 1457 !! *** ROUTINE interp1 *** 1458 1458 !! 1459 !! ** Purpose : Calculate the first order of deri avtive of1459 !! ** Purpose : Calculate the first order of derivative of 1460 1460 !! a cubic spline function y=a+b*x+c*x^2+d*x^3 1461 1461 !! -
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
r7753 r8215 37 37 38 38 ! ! Parameter to control the type of lateral viscous operator 39 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in setting the operator40 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 ! without operator (i.e. no lateral viscous trend)39 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 !: error in setting the operator 40 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 !: without operator (i.e. no lateral viscous trend) 41 41 ! !! laplacian ! bilaplacian ! 42 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 ! iso-level operator43 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 ! iso-neutral or geopotential operator42 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 !: iso-level operator 43 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 !: iso-neutral or geopotential operator 44 44 45 INTEGER :: nldf !type of lateral diffusion used defined from ln_dynldf_... (namlist logicals)45 INTEGER, PUBLIC :: nldf !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 46 46 47 47 !! * Substitutions -
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r6140 r8215 37 37 PUBLIC dyn_ldf_iso_alloc ! called by nemogcm.F90 38 38 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akzu, akzv !: vertical component of rotated lateral viscosity 40 39 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso) 40 42 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v ! - - … … 53 55 !! *** ROUTINE dyn_ldf_iso_alloc *** 54 56 !!---------------------------------------------------------------------- 55 ALLOCATE( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , &56 & zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )57 ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , & 58 & akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc ) 57 59 ! 58 60 IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.') … … 99 101 !! 100 102 !! ** Action : 101 !! Update (ua,va) arrays with the before geopotential biharmonic 102 !! mixing trend. 103 !! Update (avmu,avmv) to accompt for the diagonal vertical component 104 !! of the rotated operator in dynzdf module 103 !! -(ua,va) updated with the before geopotential harmonic mixing trend 104 !! -(akzu,akzv) to accompt for the diagonal vertical component 105 !! of the rotated operator in dynzdf module 105 106 !!---------------------------------------------------------------------- 106 107 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 144 145 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 145 146 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 146 147 !!bug 148 IF( kt == nit000 ) then 149 IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)), & 150 & ' wi', sqrt(MAXVAL(wslpi)) , ' wj', sqrt(MAXVAL(wslpj)) 151 endif 152 !!end 153 ENDIF 147 ! 148 ENDIF 154 149 155 150 ! ! =============== … … 365 360 + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) & 366 361 +zdj1u(ji,jk ) + zdju (ji ,jk ) ) 367 ! update avmu (add isopycnal vertical coefficient to avmu)368 ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0369 a vmu(ji,jj,jk) = avmu(ji,jj,jk) +( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0362 ! vertical mixing coefficient (akzu) 363 ! Note: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 364 akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0 370 365 END DO 371 366 END DO … … 391 386 & + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) & 392 387 & +zdjv (ji,jk ) + zdj1v(ji ,jk ) ) 393 ! update avmv (add isopycnal vertical coefficient to avmv)394 ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0395 a vmv(ji,jj,jk) = avmv(ji,jj,jk) +( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0388 ! vertical mixing coefficient (akzv) 389 ! Note: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 390 akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0 396 391 END DO 397 392 END DO -
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r7831 r8215 16 16 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 17 17 !! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence 18 !! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90) 18 19 !!--------------------------------------------------------------------- 19 20 … … 27 28 USE dom_oce ! ocean space and time domain 28 29 USE sbc_oce ! surface boundary condition: ocean 29 USE zdf_oce ! Bottom friction coefts 30 USE zdf_oce ! vertical physics: variables 31 USE zdfdrg ! vertical physics: top/bottom drag coef. 30 32 USE sbcisf ! ice shelf variable (fwfisf) 31 33 USE sbcapr ! surface boundary condition: atmospheric pressure … … 40 42 USE updtide ! tide potential 41 43 USE sbcwave ! surface wave 44 USE diatmb ! Top,middle,bottom output 45 #if defined key_agrif 46 USE agrif_opa_interp ! agrif 47 #endif 48 #if defined key_asminc 49 USE asminc ! Assimilation increment 50 #endif 42 51 ! 43 52 USE in_out_manager ! I/O manager … … 47 56 USE iom ! IOM library 48 57 USE restart ! only for lrst_oce 49 USE wrk_nemo ! Memory Allocation50 58 USE timing ! Timing 51 USE diatmb ! Top,middle,bottom output52 #if defined key_agrif53 USE agrif_opa_interp ! agrif54 #endif55 #if defined key_asminc56 USE asminc ! Assimilation increment57 #endif58 59 59 60 60 IMPLICIT NONE … … 66 66 PUBLIC ts_rst ! " " " " 67 67 68 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro69 REAL(wp),SAVE :: rdtbt ! Barotropic time step70 71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 !: 1st & 2nd weights used in time filtering of barotropic fields72 73 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff_f/h at F points74 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne !: triad of coriolis parameter75 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse !: (only used with een vorticity scheme)76 77 68 !! Time filtered arrays at baroclinic time step: 78 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 69 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 70 71 INTEGER , SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro 72 REAL(wp), SAVE :: rdtbt ! Barotropic time step 73 ! 74 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields 75 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points 76 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter 77 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 78 79 REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios 80 REAL(wp) :: r1_8 = 0.125_wp ! 81 REAL(wp) :: r1_4 = 0.25_wp ! 82 REAL(wp) :: r1_2 = 0.5_wp ! 79 83 80 84 !! * Substitutions 81 85 # include "vectopt_loop_substitute.h90" 82 86 !!---------------------------------------------------------------------- 83 !! NEMO/OPA 3.5 , NEMO Consortium (2013)87 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 84 88 !! $Id$ 85 89 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 137 141 INTEGER, INTENT(in) :: kt ! ocean time-step index 138 142 ! 139 LOGICAL :: ll_fw_start ! if true, forward integration140 LOGICAL :: ll_init ! if true, special startup of 2d equations141 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D142 143 INTEGER :: ji, jj, jk, jn ! dummy loop indices 143 INTEGER :: ikbu, ikbv, noffset ! local integers 144 INTEGER :: iktu, iktv ! local integers 145 REAL(wp) :: zmdi 146 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 148 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 149 REAL(wp) :: zu_spg, zv_spg ! - - 150 REAL(wp) :: zhura, zhvra ! - - 151 REAL(wp) :: za0, za1, za2, za3 ! - - 152 ! 153 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 154 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 155 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 156 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 157 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 158 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 159 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 144 LOGICAL :: ll_fw_start ! =T : forward integration 145 LOGICAL :: ll_init ! =T : special startup of 2d equations 146 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 147 INTEGER :: ikbu, iktu, noffset ! local integers 148 INTEGER :: ikbv, iktv ! - - 149 REAL(wp) :: z1_2dt_b, z2dt_bf ! local scalars 150 REAL(wp) :: zx1, zx2, zu_spg, zhura ! - - 151 REAL(wp) :: zy1, zy2, zv_spg, zhvra ! - - 152 REAL(wp) :: za0, za1, za2, za3 ! - - 153 REAL(wp) :: zmdi, zztmp ! - - 154 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 155 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 156 REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv 157 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e 158 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 159 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 160 ! 161 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 160 162 !!---------------------------------------------------------------------- 161 163 ! 162 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') 163 ! 164 ! !* Allocate temporary arrays 165 CALL wrk_alloc( jpi,jpj, zsshp2_e, zhdiv ) 166 CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd) 167 CALL wrk_alloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 168 CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 169 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 170 CALL wrk_alloc( jpi,jpj, zhf ) 171 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 164 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') 165 ! 166 IF( ln_wd ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) ) 172 167 ! 173 168 zmdi=1.e+20 ! missing data indicator for masking 174 ! !* Local constant initialization 175 z1_12 = 1._wp / 12._wp 176 z1_8 = 0.125_wp 177 z1_4 = 0.25_wp 178 z1_2 = 0.5_wp 179 zraur = 1._wp / rau0 169 ! 180 170 ! ! reciprocal of baroclinic time step 181 171 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt … … 210 200 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 211 201 ! 202 ENDIF 203 ! 204 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities) 205 DO jj = 2, jpjm1 206 DO ji = fs_2, fs_jpim1 ! vector opt. 207 zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 208 zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 209 END DO 210 END DO 211 ELSE ! bottom friction only 212 DO jj = 2, jpjm1 213 DO ji = fs_2, fs_jpim1 ! vector opt. 214 zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 215 zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 216 END DO 217 END DO 212 218 ENDIF 213 219 ! … … 263 269 !!gm 264 270 !! 265 IF ( .not.ln_sco ) THEN271 IF( .NOT.ln_sco ) THEN 266 272 267 273 !!gm agree the JC comment : this should be done in a much clear way … … 314 320 IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 315 321 ll_fw_start=.FALSE. 316 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)322 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 317 323 ENDIF 318 324 … … 363 369 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 364 370 ! energy conserving formulation for planetary vorticity term 365 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )366 zv_trd(ji,jj) = -z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )371 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 372 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 367 373 END DO 368 374 END DO … … 371 377 DO jj = 2, jpjm1 372 378 DO ji = fs_2, fs_jpim1 ! vector opt. 373 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &379 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 374 380 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 375 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &381 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 376 382 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 377 383 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) … … 383 389 DO jj = 2, jpjm1 384 390 DO ji = fs_2, fs_jpim1 ! vector opt. 385 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &391 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 386 392 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 387 393 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 388 394 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 389 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &395 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 390 396 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 391 397 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & … … 399 405 ! ! ---------------------------------------------------- 400 406 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 401 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters402 DO jj = 2, jpjm1403 DO ji = 2, jpim1404 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > &405 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &406 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &407 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 408 DO jj = 2, jpjm1 409 DO ji = 2, jpim1 410 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 411 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. & 412 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 407 413 & > rn_wdmin1 + rn_wdmin2 408 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 409 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 410 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 411 412 IF(ll_tmp1) THEN 413 zcpx(ji,jj) = 1.0_wp 414 ELSE IF(ll_tmp2) THEN 415 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 416 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 417 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 418 ELSE 419 zcpx(ji,jj) = 0._wp 420 END IF 421 422 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 414 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 415 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 416 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 417 ! 418 IF(ll_tmp1) THEN 419 zcpx(ji,jj) = 1.0_wp 420 ELSE IF(ll_tmp2) THEN ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 421 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 422 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 423 ELSE 424 zcpx(ji,jj) = 0._wp 425 ENDIF 426 ! 427 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 423 428 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. & 424 429 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 425 430 & > rn_wdmin1 + rn_wdmin2 426 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( &431 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 427 432 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 428 433 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 429 430 IF(ll_tmp1) THEN431 zcpy(ji,jj) = 1.0_wp432 ELSE IF(ll_tmp2) THEN433 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here434 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &435 &/ (sshn(ji,jj+1) - sshn(ji,jj )) )436 ELSE437 zcpy(ji,jj) = 0._wp438 ENDIF439 END DO440 END DO441 442 DO jj = 2, jpjm1443 DO ji = 2, jpim1444 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &445 &* r1_e1u(ji,jj) * zcpx(ji,jj)446 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &447 &* r1_e2v(ji,jj) * zcpy(ji,jj)448 END DO449 END DO450 434 ! 435 IF(ll_tmp1) THEN 436 zcpy(ji,jj) = 1.0_wp 437 ELSE IF(ll_tmp2) THEN 438 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 439 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 440 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 441 ELSE 442 zcpy(ji,jj) = 0._wp 443 ENDIF 444 END DO 445 END DO 446 ! 447 DO jj = 2, jpjm1 448 DO ji = 2, jpim1 449 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 450 & * r1_e1u(ji,jj) * zcpx(ji,jj) 451 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 452 & * r1_e2v(ji,jj) * zcpy(ji,jj) 453 END DO 454 END DO 455 ! 451 456 ELSE 452 453 DO jj = 2, jpjm1454 DO ji = fs_2, fs_jpim1 ! vector opt.455 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj)456 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj)457 END DO458 END DO459 ENDIF460 461 ENDIF 462 457 ! 458 DO jj = 2, jpjm1 459 DO ji = fs_2, fs_jpim1 ! vector opt. 460 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 461 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 462 END DO 463 END DO 464 ENDIF 465 ! 466 ENDIF 467 ! 463 468 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 464 469 DO ji = fs_2, fs_jpim1 … … 468 473 END DO 469 474 ! 470 ! ! Add bottomstress contribution from baroclinic velocities:471 IF (ln_bt_fw) THEN475 ! ! Add BOTTOM stress contribution from baroclinic velocities: 476 IF( ln_bt_fw ) THEN 472 477 DO jj = 2, jpjm1 473 478 DO ji = fs_2, fs_jpim1 ! vector opt. … … 491 496 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 492 497 IF( ln_wd ) THEN 493 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 494 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 498 zztmp = - 1._wp / rdtbt 499 DO jj = 2, jpjm1 500 DO ji = fs_2, fs_jpim1 ! vector opt. 501 zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) * zwx(ji,jj) 502 zv_frc(ji,jj) = zv_frc(ji,jj) + MAX( r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) * zwy(ji,jj) 503 END DO 504 END DO 495 505 ELSE 496 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 497 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 506 DO jj = 2, jpjm1 507 DO ji = fs_2, fs_jpim1 ! vector opt. 508 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj) 509 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj) 510 END DO 511 END DO 498 512 END IF 499 513 ! 500 ! ! Add top stress contribution from baroclinic velocities: 501 IF( ln_bt_fw ) THEN 514 IF( ln_isfcav ) THEN ! Add TOP stress contribution from baroclinic velocities: 515 IF( ln_bt_fw ) THEN 516 DO jj = 2, jpjm1 517 DO ji = fs_2, fs_jpim1 ! vector opt. 518 iktu = miku(ji,jj) 519 iktv = mikv(ji,jj) 520 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 521 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 522 END DO 523 END DO 524 ELSE 525 DO jj = 2, jpjm1 526 DO ji = fs_2, fs_jpim1 ! vector opt. 527 iktu = miku(ji,jj) 528 iktv = mikv(ji,jj) 529 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 530 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 531 END DO 532 END DO 533 ENDIF 534 ! 535 ! Note that the "unclipped" top friction parameter is used even with explicit drag 536 DO jj = 2, jpjm1 537 DO ji = fs_2, fs_jpim1 ! vector opt. 538 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj) 539 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj) 540 END DO 541 END DO 542 ENDIF 543 ! 544 IF( ln_bt_fw ) THEN ! Add wind forcing 502 545 DO jj = 2, jpjm1 503 546 DO ji = fs_2, fs_jpim1 ! vector opt. 504 iktu = miku(ji,jj) 505 iktv = mikv(ji,jj) 506 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 507 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 547 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj) 548 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj) 508 549 END DO 509 550 END DO 510 551 ELSE 552 zztmp = r1_rau0 * r1_2 511 553 DO jj = 2, jpjm1 512 554 DO ji = fs_2, fs_jpim1 ! vector opt. 513 iktu = miku(ji,jj) 514 iktv = mikv(ji,jj) 515 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 516 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 517 END DO 518 END DO 519 ENDIF 520 ! 521 ! Note that the "unclipped" top friction parameter is used even with explicit drag 522 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * tfrua(:,:) * zwx(:,:) 523 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * tfrva(:,:) * zwy(:,:) 524 ! 525 IF (ln_bt_fw) THEN ! Add wind forcing 526 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:) 527 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:) 528 ELSE 529 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:) 530 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:) 555 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj) 556 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj) 557 END DO 558 END DO 531 559 ENDIF 532 560 ! 533 IF ( ln_apr_dyn ) THEN! Add atm pressure forcing534 IF (ln_bt_fw) THEN561 IF( ln_apr_dyn ) THEN ! Add atm pressure forcing 562 IF( ln_bt_fw ) THEN 535 563 DO jj = 2, jpjm1 536 564 DO ji = fs_2, fs_jpim1 ! vector opt. … … 542 570 END DO 543 571 ELSE 572 zztmp = grav * r1_2 544 573 DO jj = 2, jpjm1 545 574 DO ji = fs_2, fs_jpim1 ! vector opt. 546 zu_spg = grav * z1_2* ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &547 & 548 zv_spg = grav * z1_2* ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &549 & 575 zu_spg = zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 576 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 577 zv_spg = zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 578 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 550 579 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 551 580 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 558 587 ! ! Surface net water flux and rivers 559 588 IF (ln_bt_fw) THEN 560 zssh_frc(:,:) = zraur* ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )589 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 561 590 ELSE 562 zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 563 & + fwfisf(:,:) + fwfisf_b(:,:) ) 591 zztmp = r1_rau0 * r1_2 592 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 593 & + fwfisf(:,:) + fwfisf_b(:,:) ) 564 594 ENDIF 565 595 ! … … 657 687 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 658 688 DO ji = 2, fs_jpim1 ! Vector opt. 659 zwx(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &689 zwx(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 660 690 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 661 691 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 662 zwy(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &692 zwy(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 663 693 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 664 694 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) … … 734 764 DO jj = 2, jpjm1 735 765 DO ji = 2, jpim1 ! NO Vector Opt. 736 zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &766 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 737 767 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 738 768 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 739 zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &769 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 740 770 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 741 771 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) … … 813 843 DO jj = 2, jpjm1 814 844 DO ji = 2, jpim1 815 zx1 = z1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) &845 zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 816 846 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 817 847 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 818 zy1 = z1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) &848 zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 819 849 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 820 850 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) … … 840 870 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 841 871 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 842 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )843 zv_trd(ji,jj) =- z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )872 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 873 zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 844 874 END DO 845 875 END DO … … 848 878 DO jj = 2, jpjm1 849 879 DO ji = fs_2, fs_jpim1 ! vector opt. 850 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &880 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 851 881 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 852 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &882 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 853 883 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 854 884 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) … … 860 890 DO jj = 2, jpjm1 861 891 DO ji = fs_2, fs_jpim1 ! vector opt. 862 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &892 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 863 893 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 864 894 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 865 895 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 866 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &896 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 867 897 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 868 898 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & … … 885 915 ENDIF 886 916 ! 887 ! Add bottom stresses: 888 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 889 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 890 ! 891 ! Add top stresses: 892 zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:) 893 zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 917 DO jj = 2, jpjm1 918 DO ji = fs_2, fs_jpim1 ! vector opt. 919 ! Add top/bottom stresses: 920 !!gm old/new 921 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 922 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 923 !!gm 924 END DO 925 END DO 894 926 ! 895 927 ! Surface pressure trend: … … 1025 1057 vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:) 1026 1058 ELSE 1027 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:)1028 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:)1059 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:) 1060 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:) 1029 1061 END IF 1030 1062 … … 1044 1076 DO jj = 1, jpjm1 1045 1077 DO ji = 1, jpim1 ! NO Vector Opt. 1046 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &1078 zsshu_a(ji,jj) = r1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1047 1079 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1048 1080 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1049 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &1081 zsshv_a(ji,jj) = r1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1050 1082 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1051 1083 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) … … 1091 1123 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) 1092 1124 ! 1093 CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv ) 1094 CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd ) 1095 CALL wrk_dealloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 1096 CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 1097 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 1098 CALL wrk_dealloc( jpi,jpj, zhf ) 1099 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 1125 IF( ln_wd ) DEALLOCATE( zcpx, zcpy ) 1100 1126 ! 1101 1127 IF ( ln_diatmb ) THEN … … 1248 1274 INTEGER :: ji ,jj ! dummy loop indices 1249 1275 REAL(wp) :: zxr2, zyr2, zcmax ! local scalar 1250 REAL(wp), POINTER, DIMENSION(:,:) :: zcu1276 REAL(wp), DIMENSION(jpi,jpj) :: zcu 1251 1277 !!---------------------------------------------------------------------- 1252 1278 ! 1253 1279 ! Max courant number for ext. grav. waves 1254 !1255 CALL wrk_alloc( jpi,jpj, zcu )1256 1280 ! 1257 1281 DO jj = 1, jpj … … 1320 1344 ENDIF 1321 1345 ! 1322 CALL wrk_dealloc( jpi,jpj, zcu )1323 !1324 1346 END SUBROUTINE dyn_spg_ts_init 1325 1347 -
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r7753 r8215 6 6 !! History : 1.0 ! 2005-11 (G. Madec) Original code 7 7 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 8 !! ----------------------------------------------------------------------9 10 !!---------------------------------------------------------------------- 11 !! dyn_zdf : Update the momentum trend with the vertical diffusion12 !! dyn_zdf _init : initializations of the vertical diffusion scheme8 !! 4.0 ! 2017-06 (G. Madec) remove the explicit time-stepping option + avm at t-point 9 !!---------------------------------------------------------------------- 10 11 !!---------------------------------------------------------------------- 12 !! dyn_zdf : compute the after velocity through implicit calculation of vertical mixing 13 13 !!---------------------------------------------------------------------- 14 14 USE oce ! ocean dynamics and tracers variables 15 USE phycst ! physical constants 15 16 USE dom_oce ! ocean space and time domain variables 17 USE sbc_oce ! surface boundary condition: ocean 16 18 USE zdf_oce ! ocean vertical physics variables 17 USE dynzdf_exp ! vertical diffusion: explicit (dyn_zdf_exp routine) 18 USE dynzdf_imp ! vertical diffusion: implicit (dyn_zdf_imp routine) 19 USE zdfdrg ! vertical physics: top/bottom drag coef. 20 USE dynadv ,ONLY: ln_dynadv_vec ! dynamics: advection form 21 USE dynldf ,ONLY: nldf, np_lap_i ! dynamics: type of lateral mixing 22 USE dynldf_iso,ONLY: akzu, akzv ! dynamics: vertical component of rotated lateral mixing 19 23 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 20 24 USE trd_oce ! trends: ocean variables … … 24 28 USE lib_mpp ! MPP library 25 29 USE prtctl ! Print control 26 USE wrk_nemo ! Memory Allocation27 30 USE timing ! Timing 28 31 … … 30 33 PRIVATE 31 34 32 PUBLIC dyn_zdf ! routine called by step.F90 33 PUBLIC dyn_zdf_init ! routine called by opa.F90 34 35 INTEGER :: nzdf = 0 ! type vertical diffusion algorithm used, defined from ln_zdf... namlist logicals 35 PUBLIC dyn_zdf ! routine called by step.F90 36 37 REAL(wp) :: r_vvl ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise 36 38 37 39 !! * Substitutions 38 40 # include "vectopt_loop_substitute.h90" 39 41 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 3.3 , NEMO Consortium (2010)42 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 41 43 !! $Id$ 42 44 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 43 45 !!---------------------------------------------------------------------- 44 45 46 CONTAINS 46 47 … … 49 50 !! *** ROUTINE dyn_zdf *** 50 51 !! 51 !! ** Purpose : compute the vertical ocean dynamics physics. 52 !! ** Purpose : compute the trend due to the vert. momentum diffusion 53 !! together with the Leap-Frog time stepping using an 54 !! implicit scheme. 55 !! 56 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 57 !! ua = ub + 2*dt * ua vector form or linear free surf. 58 !! ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a otherwise 59 !! - update the after velocity with the implicit vertical mixing. 60 !! This requires to solver the following system: 61 !! ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ] 62 !! with the following surface/top/bottom boundary condition: 63 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) 64 !! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 65 !! 66 !! ** Action : (ua,va) after velocity 52 67 !!--------------------------------------------------------------------- 53 !! 54 INTEGER, INTENT( in ) :: kt ! ocean time-step index 55 ! 56 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 68 INTEGER , INTENT(in) :: kt ! ocean time-step index 69 ! 70 INTEGER :: ji, jj, jk ! dummy loop indices 71 INTEGER :: iku, ikv ! local integers 72 REAL(wp) :: zzwi, ze3ua, zdt ! local scalars 73 REAL(wp) :: zzws, ze3va ! - - 74 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwd, zws ! 3D workspace 75 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - - 57 76 !!--------------------------------------------------------------------- 58 77 ! 59 78 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf') 60 79 ! 61 ! ! set time step 80 IF( kt == nit000 ) THEN !* initialization 81 IF(lwp) WRITE(numout,*) 82 IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 83 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 84 ! 85 If( ln_linssh ) THEN ; r_vvl = 0._wp ! non-linear free surface indicator 86 ELSE ; r_vvl = 1._wp 87 ENDIF 88 ENDIF 89 ! !* set time step 62 90 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdt (restart with Euler time stepping) 63 91 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! = 2 rdt (leapfrog) 64 92 ENDIF 65 93 66 IF( l_trddyn ) THEN !temporary save of ta and sa trends67 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv)94 IF( l_trddyn ) THEN !* temporary save of ta and sa trends 95 ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 68 96 ztrdu(:,:,:) = ua(:,:,:) 69 97 ztrdv(:,:,:) = va(:,:,:) 70 98 ENDIF 71 72 SELECT CASE ( nzdf ) ! compute lateral mixing trend and add it to the general trend 73 ! 74 CASE ( 0 ) ; CALL dyn_zdf_exp( kt, r2dt ) ! explicit scheme 75 CASE ( 1 ) ; CALL dyn_zdf_imp( kt, r2dt ) ! implicit scheme 76 ! 77 END SELECT 78 99 ! 100 ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in ua,va) 101 ! 102 ! ! time stepping except vertical diffusion 103 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 104 DO jk = 1, jpkm1 105 ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 106 va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 107 END DO 108 ELSE ! applied on thickness weighted velocity 109 DO jk = 1, jpkm1 110 ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) & 111 & + r2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 112 va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) & 113 & + r2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 114 END DO 115 ENDIF 116 ! ! add top/bottom friction 117 ! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom. 118 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 119 ! not lead to the effective stress seen over the whole barotropic loop. 120 ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 121 IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 122 DO jk = 1, jpkm1 ! remove barotropic velocities 123 ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 124 va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 125 END DO 126 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only 127 DO ji = fs_2, fs_jpim1 ! vector opt. 128 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 129 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 130 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 131 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 132 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 133 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 134 END DO 135 END DO 136 IF( ln_isfcav ) THEN ! Ocean cavities (ISF) 137 DO jj = 2, jpjm1 138 DO ji = fs_2, fs_jpim1 ! vector opt. 139 iku = miku(ji,jj) ! top ocean level at u- and v-points 140 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 141 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 142 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 143 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 144 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 145 END DO 146 END DO 147 END IF 148 ENDIF 149 ! 150 ! !== Vertical diffusion on u ==! 151 ! 152 ! !* Matrix construction 153 zdt = r2dt * 0.5 154 IF( nldf == np_lap_i ) THEN ! rotated lateral mixing: add its vertical mixing (akzu) 155 DO jk = 1, jpkm1 156 DO jj = 2, jpjm1 157 DO ji = fs_2, fs_jpim1 ! vector opt. 158 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 159 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 160 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 161 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 162 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 163 zwi(ji,jj,jk) = zzwi 164 zws(ji,jj,jk) = zzws 165 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 166 END DO 167 END DO 168 END DO 169 ELSE ! standard case 170 DO jk = 1, jpkm1 171 DO jj = 2, jpjm1 172 DO ji = fs_2, fs_jpim1 ! vector opt. 173 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 174 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 175 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 176 zwi(ji,jj,jk) = zzwi 177 zws(ji,jj,jk) = zzws 178 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 179 END DO 180 END DO 181 END DO 182 ENDIF 183 ! 184 DO jj = 2, jpjm1 !* Surface boundary conditions 185 DO ji = fs_2, fs_jpim1 ! vector opt. 186 zwi(ji,jj,1) = 0._wp 187 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 188 END DO 189 END DO 190 ! 191 ! !== Apply semi-implicit bottom friction ==! 192 ! 193 ! Only needed for semi-implicit bottom friction setup. The explicit 194 ! bottom friction has been included in "u(v)a" which act as the R.H.S 195 ! column vector of the tri-diagonal matrix equation 196 ! 197 IF ( ln_drgimp ) THEN ! implicit bottom friction 198 DO jj = 2, jpjm1 199 DO ji = 2, jpim1 200 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 201 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 202 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 203 END DO 204 END DO 205 IF ( ln_isfcav ) THEN ! top friction (always implicit) 206 DO jj = 2, jpjm1 207 DO ji = 2, jpim1 208 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 209 iku = miku(ji,jj) ! ocean top level at u- and v-points 210 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 211 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 212 END DO 213 END DO 214 END IF 215 ENDIF 216 ! 217 ! Matrix inversion starting from the first level 218 !----------------------------------------------------------------------- 219 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 220 ! 221 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 222 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 223 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 224 ! ( ... )( ... ) ( ... ) 225 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 226 ! 227 ! m is decomposed in the product of an upper and a lower triangular matrix 228 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 229 ! The solution (the after velocity) is in ua 230 !----------------------------------------------------------------------- 231 ! 232 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 233 DO jj = 2, jpjm1 234 DO ji = fs_2, fs_jpim1 ! vector opt. 235 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 236 END DO 237 END DO 238 END DO 239 ! 240 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 241 DO ji = fs_2, fs_jpim1 ! vector opt. 242 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 243 ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 244 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 245 END DO 246 END DO 247 DO jk = 2, jpkm1 248 DO jj = 2, jpjm1 249 DO ji = fs_2, fs_jpim1 250 ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 251 END DO 252 END DO 253 END DO 254 ! 255 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 256 DO ji = fs_2, fs_jpim1 ! vector opt. 257 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 258 END DO 259 END DO 260 DO jk = jpk-2, 1, -1 261 DO jj = 2, jpjm1 262 DO ji = fs_2, fs_jpim1 263 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 264 END DO 265 END DO 266 END DO 267 ! 268 ! !== Vertical diffusion on v ==! 269 ! 270 ! !* Matrix construction 271 zdt = r2dt * 0.5 272 IF( nldf == np_lap_i ) THEN ! rotated lateral mixing: add its vertical mixing (akzu) 273 DO jk = 1, jpkm1 274 DO jj = 2, jpjm1 275 DO ji = fs_2, fs_jpim1 ! vector opt. 276 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 277 zzwi = - zdt * ( avm(ji,jj+1,jk )+ avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 278 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 279 zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 280 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 281 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 282 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 283 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 284 END DO 285 END DO 286 END DO 287 ELSE ! standard case 288 DO jk = 1, jpkm1 289 DO jj = 2, jpjm1 290 DO ji = fs_2, fs_jpim1 ! vector opt. 291 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 292 zzwi = - zdt * ( avm(ji,jj+1,jk )+ avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 293 zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 294 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 295 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 296 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 297 END DO 298 END DO 299 END DO 300 ENDIF 301 ! 302 DO jj = 2, jpjm1 !* Surface boundary conditions 303 DO ji = fs_2, fs_jpim1 ! vector opt. 304 zwi(ji,jj,1) = 0._wp 305 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 306 END DO 307 END DO 308 ! !== Apply semi-implicit top/bottom friction ==! 309 ! 310 ! Only needed for semi-implicit bottom friction setup. The explicit 311 ! bottom friction has been included in "u(v)a" which act as the R.H.S 312 ! column vector of the tri-diagonal matrix equation 313 ! 314 IF( ln_drgimp ) THEN 315 DO jj = 2, jpjm1 316 DO ji = 2, jpim1 317 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 318 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 319 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 320 END DO 321 END DO 322 IF ( ln_isfcav ) THEN 323 DO jj = 2, jpjm1 324 DO ji = 2, jpim1 325 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 326 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 327 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 328 END DO 329 END DO 330 ENDIF 331 ENDIF 332 333 ! Matrix inversion 334 !----------------------------------------------------------------------- 335 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 336 ! 337 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 338 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 339 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 340 ! ( ... )( ... ) ( ... ) 341 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 342 ! 343 ! m is decomposed in the product of an upper and lower triangular matrix 344 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 345 ! The solution (after velocity) is in 2d array va 346 !----------------------------------------------------------------------- 347 ! 348 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 349 DO jj = 2, jpjm1 350 DO ji = fs_2, fs_jpim1 ! vector opt. 351 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 352 END DO 353 END DO 354 END DO 355 ! 356 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 357 DO ji = fs_2, fs_jpim1 ! vector opt. 358 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 359 va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 360 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 361 END DO 362 END DO 363 DO jk = 2, jpkm1 364 DO jj = 2, jpjm1 365 DO ji = fs_2, fs_jpim1 ! vector opt. 366 va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 367 END DO 368 END DO 369 END DO 370 ! 371 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 372 DO ji = fs_2, fs_jpim1 ! vector opt. 373 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 374 END DO 375 END DO 376 DO jk = jpk-2, 1, -1 377 DO jj = 2, jpjm1 378 DO ji = fs_2, fs_jpim1 379 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 380 END DO 381 END DO 382 END DO 383 ! 79 384 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 80 385 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 81 386 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 82 387 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 83 CALL wrk_dealloc( jpi, jpj, jpk,ztrdu, ztrdv )388 DEALLOCATE( ztrdu, ztrdv ) 84 389 ENDIF 85 390 ! ! print mean trends (used for debugging) … … 91 396 END SUBROUTINE dyn_zdf 92 397 93 94 SUBROUTINE dyn_zdf_init95 !!----------------------------------------------------------------------96 !! *** ROUTINE dyn_zdf_init ***97 !!98 !! ** Purpose : initializations of the vertical diffusion scheme99 !!100 !! ** Method : implicit (euler backward) scheme (default)101 !! explicit (time-splitting) scheme if ln_zdfexp=T102 !!----------------------------------------------------------------------103 USE zdftke104 USE zdfgls105 !!----------------------------------------------------------------------106 !107 ! Choice from ln_zdfexp read in namelist in zdfini108 IF( ln_zdfexp ) THEN ; nzdf = 0 ! use explicit scheme109 ELSE ; nzdf = 1 ! use implicit scheme110 ENDIF111 !112 ! Force implicit schemes113 IF( lk_zdftke .OR. lk_zdfgls ) nzdf = 1 ! TKE or GLS physics114 IF( ln_dynldf_iso ) nzdf = 1 ! iso-neutral lateral physics115 IF( ln_dynldf_hor .AND. ln_sco ) nzdf = 1 ! horizontal lateral physics in s-coordinate116 !117 IF(lwp) THEN ! Print the choice118 WRITE(numout,*)119 WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme'120 WRITE(numout,*) '~~~~~~~~~~~'121 IF( nzdf == 0 ) WRITE(numout,*) ' ===>> Explicit time-splitting scheme'122 IF( nzdf == 1 ) WRITE(numout,*) ' ===>> Implicit (euler backward) scheme'123 ENDIF124 !125 END SUBROUTINE dyn_zdf_init126 127 398 !!============================================================================== 128 399 END MODULE dynzdf
Note: See TracChangeset
for help on using the changeset viewer.
