- Timestamp:
- 2020-04-28T11:10:38+02:00 (4 years ago)
- File:
-
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90
r12529 r12822 5 5 !!====================================================================== 6 6 !! History : 3.7 ! 2014-01 (G. Madec, S. Masson) Original code, re-entrant laplacian 7 !! 4.0 ! 2020-04 (A. Nasser, G. Madec) Add symmetric mixing tensor 7 8 !!---------------------------------------------------------------------- 8 9 … … 25 26 PUBLIC dyn_ldf_lap ! called by dynldf.F90 26 27 PUBLIC dyn_ldf_blp ! called by dynldf.F90 27 28 !!anSYM 29 INTEGER, PUBLIC, PARAMETER :: np_dynldf_lap_rot = 1 ! div-rot laplacian 30 INTEGER, PUBLIC, PARAMETER :: np_dynldf_lap_sym = 2 ! symmetric laplacian (Griffies&Hallberg 2000) 31 INTEGER, PUBLIC, PARAMETER :: np_dynldf_lap_symN = 3 ! symmetric laplacian (cartesian) 32 33 INTEGER, PUBLIC, PARAMETER :: ln_dynldf_lap_typ = 1 ! choose type of laplacian (ideally from namelist) 34 !!anSYM 28 35 !! * Substitutions 29 36 # include "do_loop_substitute.h90" … … 44 51 !! ** Method : The Laplacian operator apply on horizontal velocity is 45 52 !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 53 !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 46 54 !! 47 55 !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv. 56 !! 57 !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/ 48 58 !!---------------------------------------------------------------------- 49 59 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 57 67 REAL(wp) :: zua, zva ! local scalars 58 68 REAL(wp), DIMENSION(jpi,jpj) :: zcur, zdiv 59 !!---------------------------------------------------------------------- 60 ! 69 REAL(wp), DIMENSION(jpi,jpj) :: zten, zshe ! tension (diagonal) and shearing (anti-diagonal) terms 70 !!---------------------------------------------------------------------- 71 ! 72 !!anSYM TO BE ADDED : reading of laplacian operator (ln_dynldf_lap_typ -> to be written nn_) shall be added in dyn_ldf_init 73 !! as the writing 74 !! and an integer as np_dynldf_lap for instance taken as argument by dyn_ldf_lap call in dyn_ldf 61 75 IF( kt == nit000 .AND. lwp ) THEN 62 76 WRITE(numout,*) 63 77 WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass 64 78 WRITE(numout,*) '~~~~~~~ ' 79 WRITE(numout,*) ' ln_dynldf_lap_typ = ', ln_dynldf_lap_typ 80 SELECT CASE( ln_dynldf_lap_typ ) ! print the choice of operator 81 CASE( np_dynldf_lap_rot ) ; WRITE(numout,*) ' ==>>> div-rot laplacian' 82 CASE( np_dynldf_lap_sym ) ; WRITE(numout,*) ' ==>>> symmetric laplacian (covariant form)' 83 CASE( np_dynldf_lap_symN) ; WRITE(numout,*) ' ==>>> symmetric laplacian (simple form)' 84 END SELECT 65 85 ENDIF 66 86 ! … … 69 89 ENDIF 70 90 ! 71 ! ! =============== 72 DO jk = 1, jpkm1 ! Horizontal slab 73 ! ! =============== 74 DO_2D_01_01 75 ! ! ahm * e3 * curl (computed from 1 to jpim1/jpjm1) 76 !!gm open question here : e3f at before or now ? probably now... 91 SELECT CASE( ln_dynldf_lap_typ ) 92 ! 93 CASE ( np_dynldf_lap_rot ) !== Vorticity-Divergence form ==! 94 ! 95 DO jk = 1, jpkm1 ! Horizontal slab 96 ! 97 DO_2D_01_01 98 ! ! ahm * e3 * curl (computed from 1 to jpim1/jpjm1) 99 !!gm open question here : e3f at before or now ? probably now... 77 100 !!gm note that ahmf has already been multiplied by fmask 78 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) &79 & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) &80 & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) )81 ! ! ahm * div (computed from 2 to jpi/jpj)101 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & 102 & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) & 103 & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) ) 104 ! ! ahm * div (computed from 2 to jpi/jpj) 82 105 !!gm note that ahmt has already been multiplied by tmask 83 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & 84 & * ( e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk) & 85 & + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk) ) 86 END_2D 106 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & 107 & * ( e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk) & 108 & + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk) ) 109 END_2D 110 ! 111 DO_2D_00_00 112 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * ( & 113 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & 114 & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) 115 ! 116 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * ( & 117 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm) & 118 & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) 119 END_2D 120 ! 121 END DO ! End of slab 122 ! 123 CASE ( np_dynldf_lap_sym ) !== Symmetric form ==! (Griffies&Hallberg 2000) 124 ! 125 DO jk = 1, jpkm1 ! Horizontal slab 126 ! 127 DO_2D_01_01 128 ! ! shearing stress component (F-point) NB : ahmf has already been multiplied by fmask 129 zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) & 130 & * ( e1f(ji-1,jj-1) * r1_e2f(ji-1,jj-1) & 131 & * ( pu(ji-1,jj ,jk) * r1_e1u(ji-1,jj ) - pu(ji-1,jj-1,jk) * r1_e1u(ji-1,jj-1) ) & 132 & + e2f(ji-1,jj-1) * r1_e1f(ji-1,jj-1) & 133 & * ( pv(ji ,jj-1,jk) * r1_e2v(ji ,jj-1) - pv(ji-1,jj-1,jk) * r1_e2v(ji-1,jj-1) ) ) 134 ! ! tension stress component (T-point) NB : ahmt has already been multiplied by tmask 135 zten(ji,jj) = ahmt(ji,jj,jk) & 136 & * ( e2t(ji,jj) * r1_e1t(ji,jj) & 137 & * ( pu(ji,jj,jk) * r1_e2u(ji,jj) - pu(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) & 138 & - e1t(ji,jj) * r1_e2t(ji,jj) & 139 & * ( pv(ji,jj,jk) * r1_e1v(ji,jj) - pv(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) ) 140 END_2D 141 ! 142 DO_2D_00_00 143 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & 144 & * ( ( zten(ji+1,jj ) * e2t(ji+1,jj )*e2t(ji+1,jj ) * e3t(ji+1,jj ,jk,Kmm) & 145 & - zten(ji ,jj ) * e2t(ji ,jj )*e2t(ji ,jj ) * e3t(ji ,jj ,jk,Kmm) ) * r1_e2u(ji,jj) & 146 & + ( zshe(ji ,jj ) * e1f(ji ,jj )*e1f(ji ,jj ) * e3f(ji ,jj ,jk) & 147 & - zshe(ji ,jj-1) * e1f(ji ,jj-1)*e1f(ji ,jj-1) * e3f(ji ,jj-1,jk) ) * r1_e1u(ji,jj) ) 148 ! 149 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) & 150 & * ( ( zshe(ji ,jj ) * e2f(ji ,jj )*e2f(ji ,jj ) * e3f(ji ,jj ,jk) & 151 & - zshe(ji-1,jj ) * e2f(ji-1,jj )*e2f(ji-1,jj ) * e3f(ji-1,jj ,jk) ) * r1_e2v(ji,jj) & 152 & - ( zten(ji ,jj+1) * e1t(ji ,jj+1)*e1t(ji ,jj+1) * e3t(ji ,jj+1,jk,Kmm) & 153 & - zten(ji ,jj ) * e1t(ji ,jj )*e1t(ji ,jj ) * e3t(ji ,jj ,jk,Kmm) ) * r1_e1v(ji,jj) ) 154 ! 155 END_2D 156 ! 157 END DO ! End of slab 158 ! 159 CASE ( np_dynldf_lap_symN ) !== Symmetric form ==! (naive way) 160 ! 161 DO jk = 1, jpkm1 ! Horizontal slab 162 ! 163 DO_2D_01_01 164 ! ! shearing stress component (F-point) NB : ahmf has already been multiplied by fmask 165 zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) & 166 & * ( r1_e2f(ji-1,jj-1) * ( pu(ji-1,jj ,jk) - pu(ji-1,jj-1,jk) ) & 167 & + r1_e1f(ji-1,jj-1) * ( pv(ji ,jj-1,jk) - pv(ji-1,jj-1,jk) ) ) 168 ! ! tension stress component (T-point) NB : ahmt has already been multiplied by tmask 169 zten(ji,jj) = ahmt(ji,jj,jk) & 170 & * ( r1_e1t(ji,jj) * ( pu(ji,jj,jk) - pu(ji-1,jj ,jk) ) & 171 & - r1_e2t(ji,jj) * ( pv(ji,jj,jk) - pv(ji ,jj-1,jk) ) ) 172 END_2D 173 ! 174 DO_2D_00_00 175 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & 176 & * ( zten(ji+1,jj ) * e2t(ji+1,jj ) * e3t(ji+1,jj ,jk,Kmm) & 177 & - zten(ji ,jj ) * e2t(ji ,jj ) * e3t(ji ,jj ,jk,Kmm) & 178 & + zshe(ji ,jj ) * e1f(ji ,jj ) * e3f(ji ,jj ,jk) & 179 & - zshe(ji ,jj-1) * e1f(ji ,jj-1) * e3f(ji ,jj-1,jk) ) 180 ! 181 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) & 182 & * ( zshe(ji ,jj ) * e2f(ji ,jj ) * e3f(ji ,jj ,jk) & 183 & - zshe(ji-1,jj ) * e2f(ji-1,jj ) * e3f(ji-1,jj ,jk) & 184 & - zten(ji ,jj+1) * e1t(ji ,jj+1) * e3t(ji ,jj+1,jk,Kmm) & 185 & + zten(ji ,jj ) * e1t(ji ,jj ) * e3t(ji ,jj ,jk,Kmm) ) 186 ! 187 END_2D 188 ! 189 END DO ! End of slab 190 ! 191 CASE DEFAULT ! error 192 CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for ln_dynldf_lap_typ' ) 193 END SELECT 87 194 ! 88 DO_2D_00_0089 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * ( &90 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm) &91 & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) )92 !93 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * ( &94 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm) &95 & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) )96 END_2D97 ! ! ===============98 END DO ! End of slab99 ! ! ===============100 195 ! 101 196 END SUBROUTINE dyn_ldf_lap
Note: See TracChangeset
for help on using the changeset viewer.
