Changeset 1438 for trunk/NEMO/OPA_SRC/DYN
- Timestamp:
- 2009-05-11T16:34:47+02:00 (15 years ago)
- Location:
- trunk/NEMO/OPA_SRC/DYN
- Files:
-
- 7 edited
-
dynnxt.F90 (modified) (6 diffs)
-
dynspg_exp.F90 (modified) (7 diffs)
-
dynspg_flt.F90 (modified) (14 diffs)
-
dynspg_oce.F90 (modified) (1 diff)
-
dynspg_ts.F90 (modified) (27 diffs)
-
dynvor.F90 (modified) (23 diffs)
-
wzvmod.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/dynnxt.F90
r1152 r1438 4 4 !! Ocean dynamics: time stepping 5 5 !!====================================================================== 6 6 !!====================================================================== 7 !! History : OPA ! 1987-02 (P. Andrich, D. L Hostis) Original code 8 !! ! 1990-10 (C. Levy, G. Madec) 9 !! 7.0 ! 1993-03 (M. Guyon) symetrical conditions 10 !! 8.0 ! 1997-02 (G. Madec & M. Imbard) opa, release 8.0 11 !! 8.2 ! 1997-04 (A. Weaver) Euler forward step 12 !! - ! 1997-06 (G. Madec) lateral boudary cond., lbc routine 13 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 14 !! - ! 2002-10 (C. Talandier, A-M. Treguier) Open boundary cond. 15 !! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 16 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 17 !! 3.2 ! 2009-04 (G. Madec, R.Benshila)) re-introduce the vvl option 18 !!---------------------------------------------------------------------- 19 7 20 !!---------------------------------------------------------------------- 8 21 !! dyn_nxt : update the horizontal velocity from the momentum trend 9 22 !!---------------------------------------------------------------------- 10 !! * Modules used11 23 USE oce ! ocean dynamics and tracers 12 24 USE dom_oce ! ocean space and time domain … … 29 41 PRIVATE 30 42 31 !! * Accessibility32 PUBLIC dyn_nxt ! routine called by step.F90 43 PUBLIC dyn_nxt ! routine called by step.F90 44 33 45 !! * Substitutions 34 46 # include "domzgr_substitute.h90" 35 !!---------------------------------------------------------------------- 47 !!------------------------------------------------------------------------- 48 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 49 !! $Id$ 50 !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 51 !!------------------------------------------------------------------------- 36 52 37 53 CONTAINS … … 59 75 !! - Update un,vn arrays, the now horizontal velocity 60 76 !! 61 !! History :62 !! ! 87-02 (P. Andrich, D. L Hostis) Original code63 !! ! 90-10 (C. Levy, G. Madec)64 !! ! 93-03 (M. Guyon) symetrical conditions65 !! ! 97-02 (G. Madec & M. Imbard) opa, release 8.066 !! ! 97-04 (A. Weaver) Euler forward step67 !! ! 97-06 (G. Madec) lateral boudary cond., lbc routine68 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module69 !! ! 02-10 (C. Talandier, A-M. Treguier) Open boundary cond.70 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization71 !! " ! 07-07 (D. Storkey) Calls to BDY routines.72 77 !!---------------------------------------------------------------------- 73 !! * Arguments74 78 INTEGER, INTENT( in ) :: kt ! ocean time-step index 75 76 !! * Local declarations 79 !! 77 80 INTEGER :: ji, jj, jk ! dummy loop indices 78 81 REAL(wp) :: z2dt ! temporary scalar 79 REAL(wp) :: zsshun1, zsshvn1 ! temporary scalar 80 !! Variable volume 81 REAL(wp), DIMENSION(jpi,jpj) :: & ! 2D workspace 82 zsshub, zsshun, zsshua, & 83 zsshvb, zsshvn, zsshva 84 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 85 zfse3ub, zfse3un, zfse3ua, & ! 3D workspace 86 zfse3vb, zfse3vn, zfse3va 87 !!---------------------------------------------------------------------- 88 !! OPA 9.0 , LOCEAN-IPSL (2005) 89 !! $Id$ 90 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 91 !!---------------------------------------------------------------------- 82 REAL(wp) :: zue3a , zue3n , zue3b ! temporary scalar 83 REAL(wp) :: zve3a , zve3n , zve3b ! - - 84 REAL(wp) :: ze3u_b, ze3u_n, ze3u_a ! - - 85 REAL(wp) :: ze3v_b, ze3v_n, ze3v_a ! - - 86 REAL(wp) :: zuf , zvf ! - - 87 !!---------------------------------------------------------------------- 92 88 93 89 IF( kt == nit000 ) THEN … … 102 98 103 99 !! Explicit physics with thickness weighted updates 104 IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN 105 106 ! Sea surface elevation time stepping 107 ! ----------------------------------- 108 ! 109 DO jj = 1, jpjm1 110 DO ji = 1,jpim1 111 112 ! Sea Surface Height at u-point before 113 zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 114 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) & 115 & + e1t(ji+1,jj ) * e2t(ji+1,jj ) * sshbb(ji+1,jj ) ) 116 117 ! Sea Surface Height at v-point before 118 zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 119 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) & 120 & + e1t(ji ,jj+1) * e2t(ji ,jj+1) * sshbb(ji ,jj+1) ) 121 122 ! Sea Surface Height at u-point after 123 zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 124 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) & 125 & + e1t(ji+1,jj ) * e2t(ji+1,jj ) * ssha(ji+1,jj ) ) 126 127 ! Sea Surface Height at v-point after 128 zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 129 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) & 130 & + e1t(ji ,jj+1) * e2t(ji ,jj+1) * ssha(ji ,jj+1) ) 131 132 END DO 133 END DO 134 135 ! Boundaries conditions 136 CALL lbc_lnk( zsshua, 'U', 1. ) ; CALL lbc_lnk( zsshva, 'V', 1. ) 137 CALL lbc_lnk( zsshub, 'U', 1. ) ; CALL lbc_lnk( zsshvb, 'V', 1. ) 138 139 ! Scale factors at before and after time step 140 ! ------------------------------------------- 141 CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ; CALL dom_vvl_sf( zsshua, 'U', zfse3ua ) 142 CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ; CALL dom_vvl_sf( zsshva, 'V', zfse3va ) 143 144 ! Asselin filtered scale factor at now time step 145 ! ---------------------------------------------- 146 IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN 147 CALL dom_vvl_sf_ini( 'U', zfse3un ) ; CALL dom_vvl_sf_ini( 'V', zfse3vn ) 148 ELSE 149 zsshun(:,:) = atfp * ( zsshub(:,:) + zsshua(:,:) ) + atfp1 * sshu(:,:) 150 zsshvn(:,:) = atfp * ( zsshvb(:,:) + zsshva(:,:) ) + atfp1 * sshv(:,:) 151 CALL dom_vvl_sf( zsshun, 'U', zfse3un ) ; CALL dom_vvl_sf( zsshvn, 'V', zfse3vn ) 152 ENDIF 153 154 ! Thickness weighting 155 ! ------------------- 100 101 ! Lateral boundary conditions on ( ua, va ) 102 CALL lbc_lnk( ua, 'U', -1. ) ; CALL lbc_lnk( va, 'V', -1. ) 103 104 ! Next velocity 105 ! ------------- 106 #if defined key_dynspg_flt 107 ! Leap-frog time stepping already done in dynspg_flt.F routine 108 #else 109 IF( lk_vvl ) THEN ! Varying levels 156 110 DO jk = 1, jpkm1 157 DO jj = 1, jpj 158 DO ji = 1, jpi 159 ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk) 160 va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk) 161 162 un(ji,jj,jk) = un(ji,jj,jk) * fse3u(ji,jj,jk) 163 vn(ji,jj,jk) = vn(ji,jj,jk) * fse3v(ji,jj,jk) 164 165 ub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk) 166 vb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk) 167 END DO 168 END DO 169 END DO 170 171 ENDIF 172 173 ! Lateral boundary conditions on ( ua, va ) 174 CALL lbc_lnk( ua, 'U', -1. ) 175 CALL lbc_lnk( va, 'V', -1. ) 176 177 ! ! =============== 178 DO jk = 1, jpkm1 ! Horizontal slab 179 ! ! =============== 180 ! Next velocity 181 ! ------------- 182 #if defined key_dynspg_flt 183 ! Leap-frog time stepping already done in dynspg.F routine 184 #else 185 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 186 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 187 ! Leap-frog time stepping 188 ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk) 189 va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk) 190 END DO 191 END DO 192 193 IF( lk_vvl ) THEN 194 ! Unweight velocities prior to updating open boundaries. 195 196 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 197 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 198 ua(ji,jj,jk) = ua(ji,jj,jk) / fse3u(ji,jj,jk) 199 va(ji,jj,jk) = va(ji,jj,jk) / fse3v(ji,jj,jk) 200 201 un(ji,jj,jk) = un(ji,jj,jk) / fse3u(ji,jj,jk) 202 vn(ji,jj,jk) = vn(ji,jj,jk) / fse3v(ji,jj,jk) 203 204 ub(ji,jj,jk) = ub(ji,jj,jk) / zfse3ub(ji,jj,jk) 205 vb(ji,jj,jk) = vb(ji,jj,jk) / zfse3vb(ji,jj,jk) 206 END DO 207 END DO 208 209 ENDIF 111 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 112 & + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 113 & / fse3u_a(:,:,jk) * umask(:,:,jk) 114 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 115 & + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 116 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 117 END DO 118 ELSE 119 DO jk = 1, jpkm1 120 ! Leap-frog time stepping 121 ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 122 va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 123 END DO 124 ENDIF 210 125 211 126 # if defined key_obc 212 ! ! ===============213 END DO ! End of slab214 ! ! ===============215 127 ! Update (ua,va) along open boundaries (only in the rigid-lid case) 216 128 CALL obc_dyn( kt ) 217 129 218 130 IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN 219 ! Flather boundary condition :131 ! Flather boundary condition : 220 132 ! - Update sea surface height on each open boundary 221 133 ! sshn (= after ssh) for explicit case … … 223 135 ! - Correct the barotropic velocities 224 136 CALL obc_dyn_bt( kt ) 225 226 ! Boundary conditions on sshn ( after ssh)137 ! 138 ! Boundary conditions on sshn ( after ssh) 227 139 CALL lbc_lnk( sshn, 'T', 1. ) 228 229 IF(ln_ctl) THEN ! print sum trends (used for debugging) 230 CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask) 231 ENDIF 232 233 IF ( ln_vol_cst ) CALL obc_vol( kt ) 234 235 ENDIF 236 237 ! ! =============== 238 DO jk = 1, jpkm1 ! Horizontal slab 239 ! ! =============== 140 ! 141 IF( ln_vol_cst ) CALL obc_vol( kt ) 142 ! 143 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask ) 144 ENDIF 145 240 146 # elif defined key_bdy 241 ! ! ===============242 END DO ! End of slab243 ! ! ===============244 147 ! Update (ua,va) along open boundaries (for exp or ts options). 245 IF ( lk_dynspg_exp .or. lk_dynspg_ts ) THEN246 148 IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN 149 ! 247 150 CALL bdy_dyn_frs( kt ) 248 249 IF ( ln_bdy_fla ) THEN 250 251 ua_e(:,:)=0.0 252 va_e(:,:)=0.0 253 151 ! 152 IF( ln_bdy_fla ) THEN 153 ua_e(:,:) = 0.e0 154 va_e(:,:) = 0.e0 254 155 ! Set these variables for use in bdy_dyn_fla 255 156 hu_e(:,:) = hu(:,:) 256 157 hv_e(:,:) = hv(:,:) 257 258 DO jk = 1, jpkm1 259 !! Vertically integrated momentum trends 158 DO jk = 1, jpkm1 !! Vertically integrated momentum trends 260 159 ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 261 160 va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 262 161 END DO 263 264 162 DO jk = 1 , jpkm1 265 163 ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:) 266 164 va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:) 267 165 END DO 268 269 166 CALL bdy_dta_bt( kt+1, 0) 270 167 CALL bdy_dyn_fla 271 272 168 ENDIF 273 169 ! 274 170 DO jk = 1 , jpkm1 275 171 ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:) 276 172 va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:) 277 173 END DO 278 279 ENDIF 280 281 ! ! =============== 282 DO jk = 1, jpkm1 ! Horizontal slab 283 ! ! =============== 174 ENDIF 284 175 # endif 176 285 177 # if defined key_agrif 286 ! ! ===============287 END DO ! End of slab288 ! ! ===============289 ! Update (ua,va) along open boundaries (only in the rigid-lid case)290 178 CALL Agrif_dyn( kt ) 291 ! ! ===============292 DO jk = 1, jpkm1 ! Horizontal slab293 ! ! ===============294 179 # endif 295 180 #endif 296 181 297 ! Time filter and swap of dynamics arrays 298 ! ------------------------------------------ 299 IF( neuler == 0 .AND. kt == nit000 ) THEN 300 IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels 301 ! caution: don't use (:,:) for this loop 302 ! it causes optimization problems on NEC in auto-tasking 182 ! Time filter and swap of dynamics arrays 183 ! ------------------------------------------ 184 IF( neuler == 0 .AND. kt == nit000 ) THEN 185 DO jk = 1, jpkm1 186 un(:,:,jk) = ua(:,:,jk) 187 vn(:,:,jk) = va(:,:,jk) 188 END DO 189 ELSE 190 IF( lk_vvl ) THEN ! Varying levels 191 DO jk = 1, jpkm1 ! filter applied on thickness weighted velocities 303 192 DO jj = 1, jpj 304 193 DO ji = 1, jpi 305 zsshun1 = umask(ji,jj,jk) / fse3u(ji,jj,jk) 306 zsshvn1 = vmask(ji,jj,jk) / fse3v(ji,jj,jk) 307 ub(ji,jj,jk) = un(ji,jj,jk) * zsshun1 * umask(ji,jj,jk) 308 vb(ji,jj,jk) = vn(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk) 309 zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk) 310 zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk) 311 un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 * umask(ji,jj,jk) 312 vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk) 194 ze3u_a = fse3u_a(ji,jj,jk) 195 ze3v_a = fse3v_a(ji,jj,jk) 196 ze3u_n = fse3u_n(ji,jj,jk) 197 ze3v_n = fse3v_n(ji,jj,jk) 198 ze3u_b = fse3u_b(ji,jj,jk) 199 ze3v_b = fse3v_b(ji,jj,jk) 200 ! 201 zue3a = ua(ji,jj,jk) * ze3u_a 202 zve3a = va(ji,jj,jk) * ze3v_a 203 zue3n = un(ji,jj,jk) * ze3u_n 204 zve3n = vn(ji,jj,jk) * ze3v_n 205 zue3b = ub(ji,jj,jk) * ze3u_b 206 zve3b = vb(ji,jj,jk) * ze3v_b 207 ! 208 ub(ji,jj,jk) = ( atfp * ( zue3b + zue3a ) + atfp1 * zue3n ) & 209 & / ( atfp * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 210 vb(ji,jj,jk) = ( atfp * ( zve3b + zve3a ) + atfp1 * zve3n ) & 211 & / ( atfp * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 212 un(ji,jj,jk) = ua(ji,jj,jk) * umask(ji,jj,jk) 213 vn(ji,jj,jk) = va(ji,jj,jk) * vmask(ji,jj,jk) 313 214 END DO 314 215 END DO 315 ELSE ! Fixed levels 216 END DO 217 ELSE ! Fixed levels 218 DO jk = 1, jpkm1 ! filter applied on velocities 316 219 DO jj = 1, jpj 317 220 DO ji = 1, jpi 318 ! Euler (forward) time stepping 319 ub(ji,jj,jk) = un(ji,jj,jk) 320 vb(ji,jj,jk) = vn(ji,jj,jk) 221 zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 222 zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 223 ub(ji,jj,jk) = zuf 224 vb(ji,jj,jk) = zvf 321 225 un(ji,jj,jk) = ua(ji,jj,jk) 322 226 vn(ji,jj,jk) = va(ji,jj,jk) 323 227 END DO 324 228 END DO 325 ENDIF 326 ELSE 327 IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels 328 ! caution: don't use (:,:) for this loop 329 ! it causes optimization problems on NEC in auto-tasking 330 DO jj = 1, jpj 331 DO ji = 1, jpi 332 zsshun1 = umask(ji,jj,jk) / zfse3un(ji,jj,jk) 333 zsshvn1 = vmask(ji,jj,jk) / zfse3vn(ji,jj,jk) 334 ub(ji,jj,jk) = ( atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) & 335 & + atfp1 * un(ji,jj,jk) ) * zsshun1 336 vb(ji,jj,jk) = ( atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) & 337 & + atfp1 * vn(ji,jj,jk) ) * zsshvn1 338 zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk) 339 zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk) 340 un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 341 vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 342 END DO 343 END DO 344 ELSE ! Fixed levels 345 DO jj = 1, jpj 346 DO ji = 1, jpi 347 ! Leap-frog time stepping 348 ub(ji,jj,jk) = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 349 vb(ji,jj,jk) = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 350 un(ji,jj,jk) = ua(ji,jj,jk) 351 vn(ji,jj,jk) = va(ji,jj,jk) 352 END DO 353 END DO 354 ENDIF 229 END DO 355 230 ENDIF 356 ! ! ===============357 END DO ! End of slab358 ! ! ===============359 360 IF(ln_ctl) THEN361 CALL prt_ctl(tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, &362 & tab3d_2=vn, clinfo2=' Vn: ', mask2=vmask)363 231 ENDIF 364 232 … … 367 235 #endif 368 236 237 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 238 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 239 ! 369 240 END SUBROUTINE dyn_nxt 370 241 -
trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r1146 r1438 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 !! History : 9.0 ! 2005-11 (V. Garnier, G. Madec, L. Bessieres) Original code 7 !!---------------------------------------------------------------------- 6 8 #if defined key_dynspg_exp || defined key_esopa 7 9 !!---------------------------------------------------------------------- … … 11 13 !! pressure gradient in the free surface constant 12 14 !! volume case with vector optimization 13 !! exp_rst : read/write the explicit restart fields in the ocean restart file14 15 !!---------------------------------------------------------------------- 15 !! * Modules used16 16 USE oce ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain … … 31 31 PRIVATE 32 32 33 !! * Accessibility 34 PUBLIC dyn_spg_exp ! routine called by step.F90 35 PUBLIC exp_rst ! routine called by istate.F90 33 PUBLIC dyn_spg_exp ! routine called by step.F90 36 34 37 35 !! * Substitutions … … 39 37 # include "vectopt_loop_substitute.h90" 40 38 !!---------------------------------------------------------------------- 41 !! OPA 9.0 , LOCEAN-IPSL (2005)39 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 42 40 !! $Id$ 43 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 44 42 !!---------------------------------------------------------------------- 45 46 43 47 44 CONTAINS … … 62 59 !! -1- Compute the now surface pressure gradient 63 60 !! -2- Add it to the general trend 64 !! -3- Compute the horizontal divergence of velocities65 !! - the now divergence is given by :66 !! zhdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )67 !! - integrate the horizontal divergence from the bottom68 !! to the surface69 !! - apply lateral boundary conditions on zhdivn70 !! -4- Estimate the after sea surface elevation from the kinematic71 !! surface boundary condition:72 !! zssha = sshb - 2 rdt ( zhdiv + emp )73 !! - Time filter applied on now sea surface elevation to avoid74 !! the divergence of two consecutive time-steps and swap of free75 !! surface arrays to start the next time step:76 !! sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]77 !! sshn = zssha78 !! - apply lateral boundary conditions on sshn79 61 !! 80 62 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 63 !!--------------------------------------------------------------------- 64 INTEGER, INTENT( in ) :: kt ! ocean time-step index 81 65 !! 82 !! References :83 !!84 !! History :85 !! 9.0 ! 05-11 (V. Garnier, G. Madec, L. Bessieres) Original code86 !!87 !!---------------------------------------------------------------------88 !! * Arguments89 INTEGER, INTENT( in ) :: kt ! ocean time-step index90 91 !! * Local declarations92 66 INTEGER :: ji, jj, jk ! dummy loop indices 93 REAL(wp) :: z2dt, zraur, zssha ! temporary scalars94 REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary arrays95 & zhdiv96 67 !!---------------------------------------------------------------------- 97 68 … … 104 75 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 105 76 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 106 107 CALL exp_rst( nit000, 'READ' ) ! read or initialize the following fields:108 ! ! sshb, sshn109 110 77 ENDIF 111 78 112 IF( .NOT. lk_vvl ) THEN 79 ! read or estimate sea surface height and vertically integrated velocities 80 IF( lk_obc ) CALL obc_dta_bt( kt, 0 ) 113 81 114 ! 0. Initialization115 ! -----------------116 ! read or estimate sea surface height and vertically integrated velocities117 IF( lk_obc ) CALL obc_dta_bt( kt, 0)118 z2dt = 2. * rdt ! time step: leap-frog119 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest120 zraur = 1.e0 / rauw82 ! Surface pressure gradient (now) 83 DO jj = 2, jpjm1 84 DO ji = fs_2, fs_jpim1 ! vector opt. 85 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj) 86 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 87 END DO 88 END DO 121 89 122 ! 1. Surface pressure gradient (now)123 ! ----------------------------90 ! Add the surface pressure trend to the general trend 91 DO jk = 1, jpkm1 124 92 DO jj = 2, jpjm1 125 93 DO ji = fs_2, fs_jpim1 ! vector opt. 126 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)127 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)94 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 95 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 128 96 END DO 129 97 END DO 98 END DO 130 99 131 ! 2. Add the surface pressure trend to the general trend132 ! ------------------------------------------------------133 DO jk = 1, jpkm1134 DO jj = 2, jpjm1135 DO ji = fs_2, fs_jpim1 ! vector opt.136 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)137 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)138 END DO139 END DO140 END DO141 142 ! 3. Vertical integration of horizontal divergence of velocities143 ! --------------------------------144 zhdiv(:,:) = 0.e0145 DO jk = jpkm1, 1, -1146 DO jj = 2, jpjm1147 DO ji = fs_2, fs_jpim1 ! vector opt.148 zhdiv(ji,jj) = zhdiv(ji,jj) + ( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * un(ji ,jj ,jk) &149 & - e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * un(ji-1,jj ,jk) &150 & + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * vn(ji ,jj ,jk) &151 & - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * vn(ji ,jj-1,jk) ) &152 & / ( e1t(ji,jj) * e2t(ji,jj) )153 END DO154 END DO155 END DO156 157 #if defined key_obc158 ! open boundaries (div must be zero behind the open boundary)159 ! mpp remark: The zeroing of zhdiv can probably be extended to 1->jpi/jpj for the correct row/column160 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0.e0 ! east161 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0.e0 ! west162 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north163 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0.e0 ! south164 #endif165 166 ! 4. Sea surface elevation time stepping167 ! --------------------------------------168 ! Euler (forward) time stepping, no time filter169 IF( neuler == 0 .AND. kt == nit000 ) THEN170 DO jj = 1, jpj171 DO ji = 1, jpi172 ! after free surface elevation173 zssha = sshb(ji,jj) - rdt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1)174 ! swap of arrays175 sshb(ji,jj) = sshn(ji,jj)176 sshn(ji,jj) = zssha177 END DO178 END DO179 ELSE180 ! Leap-frog time stepping and time filter181 DO jj = 1, jpj182 DO ji = 1, jpi183 ! after free surface elevation184 zssha = sshb(ji,jj) - z2dt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1)185 ! time filter and array swap186 sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj)187 sshn(ji,jj) = zssha188 END DO189 END DO190 ENDIF191 192 ELSE !! Variable volume, ssh time-stepping already done193 194 ! read or estimate sea surface height and vertically integrated velocities195 IF( lk_obc ) CALL obc_dta_bt( kt, 0 )196 197 ! Sea surface elevation swap198 ! -----------------------------199 !200 sshbb(:,:) = sshb(:,:) ! Needed for the dynamics time-stepping201 !202 IF( neuler == 0 .AND. kt == nit000 ) THEN203 ! No time filter204 sshb(:,:) = sshn(:,:)205 sshn(:,:) = ssha(:,:)206 ELSE207 ! Time filter208 sshb(:,:) = atfp * ( sshb(:,:) + ssha(:,:) ) + atfp1 * sshn(:,:)209 sshn(:,:) = ssha(:,:)210 ENDIF211 212 ENDIF213 214 ! Boundary conditions on sshn215 IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )216 217 ! write filtered free surface arrays in restart file218 ! --------------------------------------------------219 IF( lrst_oce ) CALL exp_rst( kt, 'WRITE' )220 221 IF(ln_ctl) THEN ! print sum trends (used for debugging)222 CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask)223 ENDIF224 225 100 END SUBROUTINE dyn_spg_exp 226 101 227 SUBROUTINE exp_rst( kt, cdrw )228 !!---------------------------------------------------------------------229 !! *** ROUTINE exp_rst ***230 !!231 !! ** Purpose : Read or write explicit arrays in restart file232 !!----------------------------------------------------------------------233 INTEGER , INTENT(in) :: kt ! ocean time-step234 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag235 !236 !!----------------------------------------------------------------------237 !238 IF( TRIM(cdrw) == 'READ' ) THEN239 IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN240 CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb(:,:) )241 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn(:,:) )242 IF( neuler == 0 ) sshb(:,:) = sshn(:,:)243 ELSE244 IF( nn_rstssh == 1 ) THEN245 sshb(:,:) = 0.e0246 sshn(:,:) = 0.e0247 ENDIF248 ENDIF249 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN250 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb (:,:) )251 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn (:,:) )252 ENDIF253 !254 END SUBROUTINE exp_rst255 102 #else 256 103 !!---------------------------------------------------------------------- … … 261 108 WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt 262 109 END SUBROUTINE dyn_spg_exp 263 SUBROUTINE exp_rst( kt, cdrw ) ! Empty routine264 INTEGER , INTENT(in) :: kt ! ocean time-step265 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag266 WRITE(*,*) 'exp_rst: You should not have seen this print! error?', kt, cdrw267 END SUBROUTINE exp_rst268 110 #endif 269 111 270 112 !!====================================================================== 271 113 END MODULE dynspg_exp -
trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r1200 r1438 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 !! History 8.0 ! 98-05 (G. Roullet) free surface 7 !! ! 98-10 (G. Madec, M. Imbard) release 8.2 8 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 9 !! " " ! 02-11 (C. Talandier, A-M Treguier) Open boundaries 10 !! 9.0 ! 04-08 (C. Talandier) New trends organization 11 !! " " ! 05-11 (V. Garnier) Surface pressure gradient organization 12 !! " " ! 06-07 (S. Masson) distributed restart using iom 13 !! " " ! 05-01 (J.Chanut, A.Sellar) Calls to BDY routines. 6 !! History OPA ! 1998-05 (G. Roullet) free surface 7 !! ! 1998-10 (G. Madec, M. Imbard) release 8.2 8 !! NEMO O.1 ! 2002-08 (G. Madec) F90: Free form and module 9 !! - ! 2002-11 (C. Talandier, A-M Treguier) Open boundaries 10 !! 1.0 ! 2004-08 (C. Talandier) New trends organization 11 !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization 12 !! 2.0 ! 2006-07 (S. Masson) distributed restart using iom 13 !! - ! 2006-08 (J.Chanut, A.Sellar) Calls to BDY routines. 14 !! 3.2 ! 2009-03 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 14 15 !!---------------------------------------------------------------------- 15 16 #if defined key_dynspg_flt || defined key_esopa 16 17 !!---------------------------------------------------------------------- 17 18 !! 'key_dynspg_flt' filtered free surface 18 !!----------------------------------------------------------------------19 19 !!---------------------------------------------------------------------- 20 20 !! dyn_spg_flt : update the momentum trend with the surface pressure … … 53 53 PRIVATE 54 54 55 PUBLIC dyn_spg_flt ! routine called by step.F9056 PUBLIC flt_rst ! routine called by istate.F9055 PUBLIC dyn_spg_flt ! routine called by step.F90 56 PUBLIC flt_rst ! routine called by istate.F90 57 57 58 58 !! * Substitutions … … 60 60 # include "vectopt_loop_substitute.h90" 61 61 !!---------------------------------------------------------------------- 62 !! OPA 9.0 , LOCEAN-IPSL (2005)62 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 63 63 !! $Id$ 64 64 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 80 80 !! spgv = 1/rau0 d/dy(ps) = 1/e2v dj( sshn + rnu btda ) 81 81 !! where sshn is the free surface elevation and btda is the after 82 !! of the free surface elevation 83 !! -1- compute the after sea surface elevation from the kinematic 84 !! surface boundary condition: 85 !! zssha = sshb + 2 rdt ( wn - emp ) 86 !! Time filter applied on now sea surface elevation to avoid 87 !! the divergence of two consecutive time-steps and swap of free 88 !! surface arrays to start the next time step: 89 !! sshb = sshn + atfp * [ sshb + zssha - 2 sshn ] 90 !! sshn = zssha 91 !! -2- evaluate the surface presure trend (including the addi- 82 !! time derivative of the free surface elevation 83 !! -1- evaluate the surface presure trend (including the addi- 92 84 !! tional force) in three steps: 93 85 !! a- compute the right hand side of the elliptic equation: … … 105 97 !! iterative solver 106 98 !! c- apply the solver by a call to sol... routine 107 !! - 3- compute and add the free surface pressure gradient inclu-99 !! -2- compute and add the free surface pressure gradient inclu- 108 100 !! ding the additional force used to stabilize the equation. 109 101 !! … … 112 104 !! References : Roullet and Madec 1999, JGR. 113 105 !!--------------------------------------------------------------------- 114 !! * Modules used 115 USE oce , ONLY : zub => ta, & ! ta used as workspace 116 zvb => sa ! sa " " 117 118 INTEGER, INTENT( in ) :: kt ! ocean time-step index 119 INTEGER, INTENT( out ) :: kindic ! solver convergence flag (<0 if not converge) 106 USE oce, ONLY : zub => ta ! ta used as workspace 107 USE oce, ONLY : zvb => sa ! ta used as workspace 108 !! 109 INTEGER, INTENT( in ) :: kt ! ocean time-step index 110 INTEGER, INTENT( out ) :: kindic ! solver convergence flag (<0 if not converge) 120 111 !! 121 112 INTEGER :: ji, jj, jk ! dummy loop indices 122 REAL(wp) :: z2dt, z2dtg, zraur, znugdt, & ! temporary scalars 123 & znurau, zgcb, zbtd, & ! " " 124 & ztdgu, ztdgv ! " " 125 !! Variable volume 126 REAL(wp), DIMENSION(jpi,jpj) :: & 127 zsshub, zsshua, zsshvb, zsshva, zssha ! 2D workspace 128 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & ! 3D workspace 129 zfse3ub, zfse3ua, zfse3vb, zfse3va ! 3D workspace 113 REAL(wp) :: z2dt, z2dtg, zraur, znugdt ! temporary scalars 114 REAL(wp) :: znurau, zgcb, zbtd ! " " 115 REAL(wp) :: ztdgu, ztdgv ! " " 130 116 !!---------------------------------------------------------------------- 131 117 ! … … 143 129 ! when using agrif, sshn, gcx have to be read in istate 144 130 IF (.NOT. lk_agrif) CALL flt_rst( nit000, 'READ' ) ! read or initialize the following fields: 145 ! ! gcx, gcxb , sshb, sshn131 ! ! gcx, gcxb 146 132 ENDIF 147 133 … … 158 144 IF( lk_vvl ) THEN ! variable volume 159 145 160 DO jj = 1, jpjm1161 DO ji = 1,jpim1162 163 ! Sea Surface Height at u-point before164 zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &165 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &166 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )167 168 ! Sea Surface Height at v-point before169 zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &170 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &171 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )172 173 ! Sea Surface Height at u-point after174 zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &175 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) &176 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )177 178 ! Sea Surface Height at v-point after179 zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &180 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) &181 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )182 183 END DO184 END DO185 186 ! Boundaries conditions187 CALL lbc_lnk( zsshub, 'U', 1. ) ; CALL lbc_lnk( zsshvb, 'V', 1. )188 CALL lbc_lnk( zsshua, 'U', 1. ) ; CALL lbc_lnk( zsshva, 'V', 1. )189 190 ! Scale factors at before and after time step191 ! -------------------------------------------192 CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ; CALL dom_vvl_sf( zsshua, 'U', zfse3ua )193 CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ; CALL dom_vvl_sf( zsshva, 'V', zfse3va )194 195 196 ! Thickness weighting197 ! -------------------198 DO jk = 1, jpkm1199 DO jj = 1, jpj200 DO ji = 1, jpi201 ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk)202 va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk)203 204 zub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk)205 zvb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk)206 END DO207 END DO208 END DO209 210 146 ! Evaluate the masked next velocity (effect of the additional force not included) 147 ! ------------------- (thickness weighted velocity, surface pressure gradient already included in dyn_hpg) 211 148 DO jk = 1, jpkm1 212 149 DO jj = 2, jpjm1 213 150 DO ji = fs_2, fs_jpim1 ! vector opt. 214 ua(ji,jj,jk) = ( zub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) /zfse3ua(ji,jj,jk) * umask(ji,jj,jk) 215 va(ji,jj,jk) = ( zvb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) /zfse3va(ji,jj,jk) * vmask(ji,jj,jk) 151 ua(ji,jj,jk) = ( ub(ji,jj,jk) * fse3u_b(ji,jj,jk) & 152 & + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk) ) & 153 & / fse3u_a(ji,jj,jk) * umask(ji,jj,jk) 154 va(ji,jj,jk) = ( vb(ji,jj,jk) * fse3v_b(ji,jj,jk) & 155 & + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk) ) & 156 & / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk) 216 157 END DO 217 158 END DO … … 227 168 END DO 228 169 END DO 229 230 ! Add the surface pressure trend to the general trend 170 ! 171 ! add the surface pressure trend to the general trend and 172 ! evaluate the masked next velocity (effect of the additional force not included) 231 173 DO jk = 1, jpkm1 232 174 DO jj = 2, jpjm1 233 175 DO ji = fs_2, fs_jpim1 ! vector opt. 234 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)235 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)176 ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) ) ) * umask(ji,jj,jk) 177 va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) ) ) * vmask(ji,jj,jk) 236 178 END DO 237 179 END DO 238 180 END DO 239 240 ! Evaluate the masked next velocity (effect of the additional force not included) 241 DO jk = 1, jpkm1 242 DO jj = 2, jpjm1 243 DO ji = fs_2, fs_jpim1 ! vector opt. 244 ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk) 245 va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk) 246 END DO 247 END DO 248 END DO 249 181 ! 250 182 ENDIF 251 183 … … 311 243 CALL lbc_lnk( spgv, 'V', -1. ) 312 244 313 IF( lk_vvl ) CALL sol_mat( kt ) 245 IF( lk_vvl ) CALL sol_mat( kt ) ! build the matrix at kt (vvl case only) 314 246 315 247 ! Right hand side of the elliptic equation and first guess … … 347 279 ! Relative precision (computation on one processor) 348 280 ! ------------------ 349 rnorme =0. 281 rnorme =0.e0 350 282 rnorme = SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) ) 351 283 IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain … … 413 345 ENDIF 414 346 #endif 415 ! 7.Add the trends multiplied by z2dt to the after velocity416 ! ------------------------------------------------------- ----347 ! Add the trends multiplied by z2dt to the after velocity 348 ! ------------------------------------------------------- 417 349 ! ( c a u t i o n : (ua,va) here are the after velocity not the 418 350 ! trend, the leap-frog time stepping will not … … 421 353 DO jj = 2, jpjm1 422 354 DO ji = fs_2, fs_jpim1 ! vector opt. 423 ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj)) * umask(ji,jj,jk)424 va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj)) * vmask(ji,jj,jk)355 ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk) 356 va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk) 425 357 END DO 426 358 END DO 427 359 END DO 428 429 ! Sea surface elevation time stepping430 ! -----------------------------------431 IF( lk_vvl ) THEN ! after free surface elevation432 zssha(:,:) = ssha(:,:)433 ELSE434 zssha(:,:) = sshb(:,:) + z2dt * ( wn(:,:,1) - emp(:,:) * zraur ) * tmask(:,:,1)435 ENDIF436 #if defined key_obc437 # if defined key_agrif438 IF ( Agrif_Root() ) THEN439 # endif440 zssha(:,:)=zssha(:,:)*obctmsk(:,:)441 CALL lbc_lnk(zssha,'T',1.) ! absolutly compulsory !! (jmm)442 # if defined key_agrif443 ENDIF444 # endif445 #endif446 447 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler (forward) time stepping, no time filter448 ! swap of arrays449 sshb(:,:) = sshn (:,:)450 sshn(:,:) = zssha(:,:)451 ELSE ! Leap-frog time stepping and time filter452 ! time filter and array swap453 sshb(:,:) = atfp * ( sshb(:,:) + zssha(:,:) ) + atfp1 * sshn(:,:)454 sshn(:,:) = zssha(:,:)455 ENDIF456 360 457 361 ! write filtered free surface arrays in restart file 458 362 ! -------------------------------------------------- 459 363 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 460 461 ! print sum trends (used for debugging)462 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg - ssh: ', mask1=tmask )463 364 ! 464 365 END SUBROUTINE dyn_spg_flt … … 481 382 CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) ) 482 383 CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) ) 483 CALL iom_get( numror, jpdom_autoglo, 'sshb', sshb(:,:) ) 484 CALL iom_get( numror, jpdom_autoglo, 'sshn', sshn(:,:) ) 485 IF( neuler == 0 ) THEN 486 sshb(:,:) = sshn(:,:) 487 gcxb(:,:) = gcx (:,:) 488 ENDIF 384 IF( neuler == 0 ) gcxb(:,:) = gcx (:,:) 489 385 ELSE 490 386 gcx (:,:) = 0.e0 491 387 gcxb(:,:) = 0.e0 492 IF( nn_rstssh == 1 ) THEN493 sshb(:,:) = 0.e0494 sshn(:,:) = 0.e0495 ENDIF496 388 ENDIF 497 389 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 498 390 ! Caution : extra-hallow 499 391 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 500 CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )392 CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) ) 501 393 CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) ) 502 CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:) )503 CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:) )504 394 ENDIF 505 395 ! -
trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90
r1152 r1438 43 43 !! ------------------------ 44 44 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables averaged over the barotropic loop 45 sshn_b, sshb_b, & ! sea surface heigth (now, before)46 u n_b , vn_b ! vertically integrated horizontal velocities (now)45 un_e , vn_e, & ! vertically integrated horizontal velocities (now) 46 ua_e , va_e ! vertically integrated horizontal velocities (after) 47 47 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables of the explicit barotropic loop 48 sshb_e, sshn_e, ssha_e, & ! sea surface heigth (before,now,after) 49 ua_e , va_e, & ! vertically integrated horizontal velocities (after) 50 hu_e , hv_e ! depth arrays for the barotropic solution 48 sshn_e, ssha_e, & ! sea surface heigth (now, after) 49 hu_e , hv_e ! depth arrays for the barotropic solution 51 50 #endif 52 51 -
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r1241 r1438 1 1 MODULE dynspg_ts 2 2 !!====================================================================== 3 !! *** MODULE dynspg_ts *** 4 !! Ocean dynamics: surface pressure gradient trend 5 !!====================================================================== 6 !! History : 9.0 ! 04-12 (L. Bessieres, G. Madec) Original code 7 !! " " ! 05-11 (V. Garnier, G. Madec) optimization 8 !! 9.0 ! 06-08 (S. Masson) distributed restart using iom 9 !! " ! 08-01 (R. Benshila) change averaging method 10 !! " ! 07-07 (D. Storkey) calls to BDY routines 11 !!--------------------------------------------------------------------- 3 !! History : 1.0 ! 04-12 (L. Bessieres, G. Madec) Original code 4 !! - ! 05-11 (V. Garnier, G. Madec) optimization 5 !! - ! 06-08 (S. Masson) distributed restart using iom 6 !! 2.0 ! 07-07 (D. Storkey) calls to BDY routines 7 !! - ! 08-01 (R. Benshila) change averaging method 8 !!--------------------------------------------------------------------- 12 9 #if defined key_dynspg_ts || defined key_esopa 13 10 !!---------------------------------------------------------------------- … … 19 16 !! ts_rst : read/write the time-splitting restart fields in the ocean restart file 20 17 !!---------------------------------------------------------------------- 21 !! * Modules used22 18 USE oce ! ocean dynamics and tracers 23 19 USE dom_oce ! ocean space and time domain … … 49 45 PUBLIC ts_rst ! routine called by istate.F90 50 46 51 REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne, & ! triad of coriolis parameter 52 & ftsw, ftse ! (only used with een vorticity scheme) 53 47 REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne ! triad of coriolis parameter 48 REAL(wp), DIMENSION(jpi,jpj) :: ftsw, ftse ! (only used with een vorticity scheme) 54 49 55 50 !! * Substitutions 56 51 # include "domzgr_substitute.h90" 57 52 # include "vectopt_loop_substitute.h90" 58 !!---------------------------------------------------------------------- 59 !! OPA 9.0 , LOCEAN-IPSL (2005)53 !!------------------------------------------------------------------------- 54 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 60 55 !! $Id$ 61 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt62 !!---------------------------------------------------------------------- 56 !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 57 !!------------------------------------------------------------------------- 63 58 64 59 CONTAINS … … 85 80 !! (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b 86 81 !! and zvX_b (X specifying after, now or before). 87 !! -3- Update of sea surface height from time averaged barotropic 88 !! variables. 89 !! - apply lateral boundary conditions on sshn. 90 !! -4- The new general trend becomes : 91 !! ua = ua - sum_k(ua)/H + ( zua_b - sum_k(ub) )/H 82 !! -3- The new general trend becomes : 83 !! ua = ua - sum_k(ua)/H + ( zua_e - sum_k(ub) )/H 92 84 !! 93 85 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend … … 96 88 !!--------------------------------------------------------------------- 97 89 INTEGER, INTENT( in ) :: kt ! ocean time-step index 98 99 !! * Local declarations 90 !! 100 91 INTEGER :: ji, jj, jk, jit ! dummy loop indices 101 92 INTEGER :: icycle ! temporary scalar … … 107 98 zcu, zcv, zwx, zwy, zhdiv, & ! temporary arrays 108 99 zua, zva, zub, zvb, & ! " " 109 z ssha_b, zua_b, zva_b, & ! " "110 zub_e, zvb_e, 111 zu n_e, zvn_e ! " "100 zua_e, zva_e, zssha_e, & ! " " 101 zub_e, zvb_e, zsshb_e, & ! " " 102 zu_sum, zv_sum, zssh_sum 112 103 !! Variable volume 113 104 REAL(wp), DIMENSION(jpi,jpj) :: & 114 105 zspgu_1, zspgv_1, zsshun_e, zsshvn_e ! 2D workspace 115 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfse3un_e, zfse3vn_e ! 3D workspace116 106 !!---------------------------------------------------------------------- 117 107 118 108 ! Arrays initialization 119 109 ! --------------------- 120 zua_b(:,:) = 0.e0 ; zub_e(:,:) = 0.e0 ; zun_e(:,:) = 0.e0121 zva_b(:,:) = 0.e0 ; zvb_e(:,:) = 0.e0 ; zvn_e(:,:) = 0.e0122 110 zhdiv(:,:) = 0.e0 123 111 … … 133 121 ! ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 134 122 135 ssha_e(:,:) = sshn(:,:) 136 ua_e(:,:) = un_b(:,:) 137 va_e(:,:) = vn_b(:,:) 138 hu_e(:,:) = hu(:,:) 139 hv_e(:,:) = hv(:,:) 140 123 zssha_e(:,:) = sshn(:,:) 124 zua_e (:,:) = un_e(:,:) 125 zva_e (:,:) = vn_e(:,:) 126 hu_e (:,:) = hu (:,:) 127 hv_e (:,:) = hv (:,:) 141 128 IF( ln_dynvor_een ) THEN 142 129 ftne(1,:) = 0.e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0 … … 170 157 zspgv_1(ji,jj) = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 171 158 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e2v(ji,jj) 172 END DO 159 END DO 173 160 END DO 174 161 … … 193 180 zfact2 = 0.5 * 0.5 194 181 zraur = 1. / rauw ! 1 / volumic mass of pure water 195 182 196 183 ! ----------------------------------------------------------------------------- 197 184 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) … … 215 202 zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk) 216 203 ! ! Vertically integrated transports (before) 217 zub(ji,1) = zub(ji,1) + fse3u (ji,1,jk) * ub(ji,1,jk)218 zvb(ji,1) = zvb(ji,1) + fse3v (ji,1,jk) * vb(ji,1,jk)204 zub(ji,1) = zub(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 205 zvb(ji,1) = zvb(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 219 206 ! ! Planetary vorticity transport fluxes (now) 220 207 zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk) … … 228 215 zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 229 216 ! ! Vertically integrated transports (before) 230 zub(:,:) = zub(:,:) + fse3u (:,:,jk) * ub(:,:,jk)231 zvb(:,:) = zvb(:,:) + fse3v (:,:,jk) * vb(:,:,jk)217 zub(:,:) = zub(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 218 zvb(:,:) = zvb(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 232 219 ! ! Planetary vorticity (now) 233 220 zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) … … 304 291 !---------------- 305 292 ! Number of iteration of the barotropic loop 306 icycle = 3 * nn_baro / 2293 icycle = 2 * nn_baro + 1 307 294 308 295 ! variables for the barotropic equations 309 sshb_e (:,:) = sshn_b(:,:) ! (barotropic) sea surface height (before and now) 310 sshn_e (:,:) = sshn_b(:,:) 311 zub_e (:,:) = un_b (:,:) ! barotropic transports issued from the barotropic equations (before and now) 312 zvb_e (:,:) = vn_b (:,:) 313 zun_e (:,:) = un_b (:,:) 314 zvn_e (:,:) = vn_b (:,:) 315 zssha_b(:,:) = 0.e0 316 zua_b (:,:) = 0.e0 317 zva_b (:,:) = 0.e0 318 hu_e (:,:) = hu (:,:) ! (barotropic) ocean depth at u-point 319 hv_e (:,:) = hv (:,:) ! (barotropic) ocean depth at v-point 320 IF( lk_vvl ) THEN 321 zsshun_e(:,:) = sshu (:,:) ! (barotropic) sea surface height (now) at u-point 322 zsshvn_e(:,:) = sshv (:,:) ! (barotropic) sea surface height (now) at v-point 323 zfse3un_e(:,:,:) = fse3u(:,:,:) ! (barotropic) scale factors at u-point 324 zfse3un_e(:,:,:) = fse3v(:,:,:) ! (barotropic) scale factors at v-point 296 zu_sum (:,:) = 0.e0 297 zv_sum (:,:) = 0.e0 298 zssh_sum(:,:) = 0.e0 299 hu_e (:,:) = hu (:,:) ! (barotropic) ocean depth at u-point 300 hv_e (:,:) = hv (:,:) ! (barotropic) ocean depth at v-point 301 zsshb_e (:,:) = sshn_e(:,:) ! (barotropic) sea surface height (before and now) 302 ! vertical sum 303 un_e (:,:) = 0.e0 304 vn_e (:,:) = 0.e0 305 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 306 DO jk = 1, jpkm1 307 DO ji = 1, jpij 308 un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 309 vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 310 END DO 311 END DO 312 ELSE ! No vector opt. 313 DO jk = 1, jpkm1 314 un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 315 vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 316 END DO 325 317 ENDIF 318 zub_e (:,:) = un_e(:,:) 319 zvb_e (:,:) = vn_e(:,:) 320 326 321 327 322 ! set ssh corrections to 0 … … 352 347 DO jj = 2, jpjm1 353 348 DO ji = fs_2, fs_jpim1 ! vector opt. 354 zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun_e(ji ,jj) &355 & -e2u(ji-1,jj ) * zun_e(ji-1,jj) &356 & +e1v(ji ,jj ) * zvn_e(ji ,jj) &357 & -e1v(ji ,jj-1) * zvn_e(ji ,jj-1) ) &349 zhdiv(ji,jj) = ( e2u(ji ,jj ) * un_e(ji ,jj) & 350 & -e2u(ji-1,jj ) * un_e(ji-1,jj) & 351 & +e1v(ji ,jj ) * vn_e(ji ,jj) & 352 & -e1v(ji ,jj-1) * vn_e(ji ,jj-1) ) & 358 353 & / (e1t(ji,jj)*e2t(ji,jj)) 359 354 END DO … … 370 365 371 366 #if defined key_bdy 372 373 374 375 376 367 DO jj = 1, jpj 368 DO ji = 1, jpi 369 zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 370 END DO 371 END DO 377 372 #endif 378 373 … … 381 376 DO jj = 2, jpjm1 382 377 DO ji = fs_2, fs_jpim1 ! vector opt. 383 ssha_e(ji,jj) = (sshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) &384 & + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)378 zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) & 379 & + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 385 380 END DO 386 381 END DO … … 388 383 ! evolution of the barotropic transport ( following the vorticity scheme used) 389 384 ! ---------------------------------------------------------------------------- 390 zwx(:,:) = e2u(:,:) * zun_e(:,:)391 zwy(:,:) = e1v(:,:) * zvn_e(:,:)385 zwx(:,:) = e2u(:,:) * un_e(:,:) 386 zwy(:,:) = e1v(:,:) * vn_e(:,:) 392 387 393 388 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme … … 396 391 ! surface pressure gradient 397 392 IF( lk_vvl) THEN 398 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj) &399 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)400 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &401 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)393 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 394 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 395 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 396 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 402 397 ELSE 403 398 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) … … 412 407 zcvbt =-zfact2 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) 413 408 ! after transports 414 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)415 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)409 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 410 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 416 411 END DO 417 412 END DO … … 422 417 ! surface pressure gradient 423 418 IF( lk_vvl) THEN 424 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj) &425 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)426 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &427 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)419 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 420 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 421 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 422 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 428 423 ELSE 429 424 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) … … 432 427 ! enstrophy conserving formulation for planetary vorticity term 433 428 zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 434 429 + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 435 430 zx1 =-zfact1 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 436 431 + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 437 432 zcubt = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) 438 433 zcvbt = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) 439 434 ! after transports 440 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)441 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)435 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 436 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 442 437 END DO 443 438 END DO … … 449 444 ! surface pressure gradient 450 445 IF( lk_vvl) THEN 451 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj) &452 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)453 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &454 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)446 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 447 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 448 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 449 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 455 450 ELSE 456 451 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) … … 463 458 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) 464 459 ! after transports 465 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)466 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)460 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 461 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 467 462 END DO 468 463 END DO … … 477 472 478 473 ! ... Boundary conditions on ua_e, va_e, ssha_e 479 CALL lbc_lnk( ua_e , 'U', -1. )480 CALL lbc_lnk( va_e , 'V', -1. )481 CALL lbc_lnk( ssha_e, 'T', 1. )474 CALL lbc_lnk( zua_e , 'U', -1. ) 475 CALL lbc_lnk( zva_e , 'V', -1. ) 476 CALL lbc_lnk( zssha_e, 'T', 1. ) 482 477 483 478 ! temporal sum 484 479 !------------- 485 IF( jit >= nn_baro / 2 ) THEN 486 zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:) 487 zua_b (:,:) = zua_b (:,:) + ua_e (:,:) 488 zva_b (:,:) = zva_b (:,:) + va_e (:,:) 489 ENDIF 480 zu_sum (:,:) = zu_sum (:,:) + zua_e (:,:) 481 zv_sum (:,:) = zv_sum (:,:) + zva_e (:,:) 482 zssh_sum(:,:) = zssh_sum(:,:) + zssha_e(:,:) 490 483 491 484 ! Time filter and swap of dynamics arrays 492 485 ! --------------------------------------- 493 486 IF( jit == 1 ) THEN ! Euler (forward) time stepping 494 sshb_e (:,:) = sshn_e(:,:)495 zub_e (:,:) = zun_e(:,:)496 zvb_e (:,:) = zvn_e(:,:)497 sshn_e (:,:) = ssha_e(:,:)498 zun_e (:,:) =ua_e (:,:)499 zvn_e (:,:) =va_e (:,:)487 zsshb_e(:,:) = sshn_e (:,:) 488 zub_e (:,:) = un_e (:,:) 489 zvb_e (:,:) = vn_e (:,:) 490 sshn_e (:,:) = zssha_e(:,:) 491 un_e (:,:) = zua_e (:,:) 492 vn_e (:,:) = zva_e (:,:) 500 493 ELSE ! Asselin filtering 501 sshb_e (:,:) = atfp * ( sshb_e(:,:) +ssha_e(:,:) ) + atfp1 * sshn_e(:,:)502 zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:)503 zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:)504 sshn_e (:,:) = ssha_e(:,:)505 zun_e (:,:) =ua_e (:,:)506 zvn_e (:,:) =va_e (:,:)494 zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * sshn_e(:,:) 495 zub_e (:,:) = atfp * ( zub_e (:,:) + zua_e (:,:) ) + atfp1 * un_e (:,:) 496 zvb_e (:,:) = atfp * ( zvb_e (:,:) + zva_e (:,:) ) + atfp1 * vn_e (:,:) 497 sshn_e (:,:) = zssha_e(:,:) 498 un_e (:,:) = zua_e (:,:) 499 vn_e (:,:) = zva_e (:,:) 507 500 ENDIF 508 501 509 IF( lk_vvl ) THEN ! Variable volume 502 IF( lk_vvl ) THEN ! Variable volume ! needed for BDY ??? 510 503 511 504 ! Sea surface elevation … … 528 521 529 522 ! Boundaries conditions 530 CALL lbc_lnk( zsshun_e, 'U', 1. ) ; 531 532 ! Scale factors at before and after time step523 CALL lbc_lnk( zsshun_e, 'U', 1. ) ; CALL lbc_lnk( zsshvn_e, 'V', 1. ) 524 525 ! Ocean depth at U- and V-points 533 526 ! ------------------------------------------- 534 CALL dom_vvl_sf( zsshun_e, 'U', zfse3un_e ) ; CALL dom_vvl_sf( zsshvn_e, 'V', zfse3vn_e ) 535 536 ! Ocean depth at U- and V-points 537 hu_e(:,:) = 0.e0 538 hv_e(:,:) = 0.e0 539 540 DO jk = 1, jpk 541 hu_e(:,:) = hu_e(:,:) + zfse3un_e(:,:,jk) * umask(:,:,jk) 542 hv_e(:,:) = hv_e(:,:) + zfse3vn_e(:,:,jk) * vmask(:,:,jk) 543 END DO 544 545 ENDIF ! End variable volume case 546 527 hu_e(:,:) = hu_0(:,:) + zsshun_e(:,:) 528 hv_e(:,:) = hv_0(:,:) + zsshvn_e(:,:) 529 530 ! 531 ENDIF 547 532 ! ! ==================== ! 548 533 END DO ! end loop ! 549 534 ! ! ==================== ! 550 535 551 552 536 ! Time average of after barotropic variables 553 zcoef = 1.e0 / ( nn_baro + 1 ) 554 zssha_b(:,:) = zcoef * zssha_b(:,:) 555 zua_b (:,:) = zcoef * zua_b (:,:) 556 zva_b (:,:) = zcoef * zva_b (:,:) 537 zcoef = 1.e0 / ( 2 * nn_baro + 1 ) 538 un_e (:,:) = zcoef * zu_sum (:,:) 539 vn_e (:,:) = zcoef * zv_sum (:,:) 540 sshn_e(:,:) = zcoef * zssh_sum(:,:) 541 557 542 #if defined key_obc 558 543 IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) … … 561 546 IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:) 562 547 #endif 563 564 565 ! ---------------------------------------------------------------------------566 ! Phase 3 : Update sea surface height from time averaged barotropic variables567 ! ---------------------------------------------------------------------------568 569 570 ! Horizontal divergence of time averaged barotropic transports571 !-------------------------------------------------------------572 IF( .NOT. lk_vvl ) THEN573 DO jj = 2, jpjm1574 DO ji = fs_2, fs_jpim1 ! vector opt.575 zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj ) * un_b(ji-1,jj ) &576 & +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji ,jj-1) * vn_b(ji ,jj-1) ) &577 & / ( e1t(ji,jj) * e2t(ji,jj) )578 END DO579 END DO580 ENDIF581 582 #if defined key_obc && ! defined key_vvl583 ! open boundaries (div must be zero behind the open boundary)584 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column585 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east586 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west587 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north588 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south589 #endif590 591 #if defined key_bdy592 DO jj = 1, jpj593 DO ji = 1, jpi594 zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)595 END DO596 END DO597 #endif598 599 ! sea surface height600 !-------------------601 IF( lk_vvl ) THEN602 sshbb(:,:) = sshb(:,:)603 sshb (:,:) = sshn(:,:)604 sshn (:,:) = ssha(:,:)605 ELSE606 sshb (:,:) = sshn(:,:)607 sshn(:,:) = ( sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)608 ENDIF609 610 ! ... Boundary conditions on sshn611 IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )612 613 548 614 549 ! ----------------------------------------------------------------------------- 615 ! Phase 4. Coupling between general trend and barotropic estimates - (2nd step)550 ! Phase 3. Coupling between general trend and barotropic estimates - (2nd step) 616 551 ! ----------------------------------------------------------------------------- 617 552 618 ! Swap on time averaged barotropic variables 619 ! ------------------------------------------ 620 sshb_b(:,:) = sshn_b (:,:) 621 IF ( neuler == 0 .AND. kt == nit000 ) zssha_b(:,:) = sshn(:,:) 622 sshn_b(:,:) = zssha_b(:,:) 623 un_b (:,:) = zua_b (:,:) 624 vn_b (:,:) = zva_b (:,:) 625 553 554 626 555 ! add time averaged barotropic coriolis and surface pressure gradient 627 556 ! terms to the general momentum trend 628 557 ! -------------------------------------------------------------------- 629 558 DO jk=1,jpkm1 630 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( zua_b(:,:) - zub(:,:) ) / z2dt_b631 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( zva_b(:,:) - zvb(:,:) ) / z2dt_b559 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( un_e(:,:) - zub(:,:) ) / z2dt_b 560 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( vn_e(:,:) - zvb(:,:) ) / z2dt_b 632 561 END DO 633 562 … … 636 565 IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) 637 566 638 ! print sum trends (used for debugging)639 IF( ln_ctl ) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask )640 567 ! 641 568 END SUBROUTINE dyn_spg_ts … … 655 582 ! 656 583 IF( TRIM(cdrw) == 'READ' ) THEN 657 IF( iom_varid( numror, 'sshn ', ldstop = .FALSE. ) > 0 ) THEN658 CALL iom_get( numror, jpdom_autoglo, 'ssh b' , sshb(:,:) )659 CALL iom_get( numror, jpdom_autoglo, ' sshn' , sshn(:,:) )660 IF( neuler == 0 ) sshb(:,:) = sshn(:,:)584 IF( iom_varid( numror, 'sshn_e', ldstop = .FALSE. ) > 0 ) THEN 585 CALL iom_get( numror, jpdom_autoglo, 'sshn_e', sshn_e(:,:) ) ! free surface and 586 CALL iom_get( numror, jpdom_autoglo, 'un_e' , un_e (:,:) ) ! horizontal transports issued 587 CALL iom_get( numror, jpdom_autoglo, 'vn_e' , vn_e (:,:) ) ! from barotropic loop 661 588 ELSE 662 IF( nn_rstssh == 1 ) THEN 663 sshb(:,:) = 0.e0 664 sshn(:,:) = 0.e0 665 ENDIF 666 ENDIF 667 IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN 668 CALL iom_get( numror, jpdom_autoglo, 'sshb_b', sshb_b(:,:) ) ! free surface issued 669 CALL iom_get( numror, jpdom_autoglo, 'sshn_b', sshn_b(:,:) ) ! from time-splitting loop 670 CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! horizontal transports issued 671 CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 672 IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:) 673 ELSE 674 sshb_b(:,:) = sshb(:,:) 675 sshn_b(:,:) = sshn(:,:) 676 un_b (:,:) = 0.e0 677 vn_b (:,:) = 0.e0 589 sshn_e(:,:) = sshn(:,:) 590 un_e (:,:) = 0.e0 591 vn_e (:,:) = 0.e0 678 592 ! vertical sum 679 593 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 680 594 DO jk = 1, jpkm1 681 595 DO ji = 1, jpij 682 un_ b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)683 vn_ b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)596 un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 597 vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 684 598 END DO 685 599 END DO 686 600 ELSE ! No vector opt. 687 601 DO jk = 1, jpkm1 688 un_ b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)689 vn_ b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)602 un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 603 vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 690 604 END DO 691 605 ENDIF 692 606 ENDIF 693 607 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 694 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb (:,:) ) 695 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn (:,:) ) 696 CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) ) ! free surface issued 697 CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) ) ! from barotropic loop 698 CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! horizontal transports issued 699 CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 608 CALL iom_rstput( kt, nitrst, numrow, 'sshn_e', sshn_e(:,:) ) ! free surface and 609 CALL iom_rstput( kt, nitrst, numrow, 'un_e' , un_e (:,:) ) ! horizontal transports issued 610 CALL iom_rstput( kt, nitrst, numrow, 'vn_e' , vn_e (:,:) ) ! from barotropic loop 700 611 ENDIF 701 612 ! -
trunk/NEMO/OPA_SRC/DYN/dynvor.F90
r1152 r1438 5 5 !! planetary vorticity trends 6 6 !!====================================================================== 7 !! History : 1.0 ! 89-12 (P. Andrich) vor_ens: Original code 8 !! 5.0 ! 91-11 (G. Madec) vor_ene, vor_mix: Original code 9 !! 6.0 ! 96-01 (G. Madec) s-coord, suppress work arrays 10 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 11 !! 8.5 ! 04-02 (G. Madec) vor_een: Original code 12 !! 9.0 ! 03-08 (G. Madec) vor_ctl: Original code 13 !! 9.0 ! 05-11 (G. Madec) dyn_vor: Original code (new step architecture) 14 !! 9.0 ! 06-11 (G. Madec) flux form advection: add metric term 7 !! History : OPA ! 1989-12 (P. Andrich) vor_ens: Original code 8 !! 5.0 ! 1991-11 (G. Madec) vor_ene, vor_mix: Original code 9 !! 6.0 ! 1996-01 (G. Madec) s-coord, suppress work arrays 10 !! 8.5 ! 2002-08 (G. Madec) F90: Free form and module 11 !! NEMO 1.0 ! 2004-02 (G. Madec) vor_een: Original code 12 !! - ! 2003-08 (G. Madec) add vor_ctl 13 !! - ! 2005-11 (G. Madec) add dyn_vor (new step architecture) 14 !! 2.0 ! 2006-11 (G. Madec) flux form advection: add metric term 15 !! 3.2 ! 2009-04 (R. Benshila) vvl: correction of een scheme 15 16 !!---------------------------------------------------------------------- 16 17 … … 37 38 PUBLIC dyn_vor ! routine called by step.F90 38 39 39 ! !* Namelist nam_dynvor: vorticity term40 ! !!* Namelist nam_dynvor: vorticity term 40 41 LOGICAL, PUBLIC :: ln_dynvor_ene = .FALSE. !: energy conserving scheme 41 42 LOGICAL, PUBLIC :: ln_dynvor_ens = .TRUE. !: enstrophy conserving scheme … … 52 53 # include "vectopt_loop_substitute.h90" 53 54 !!---------------------------------------------------------------------- 54 !! OPA 9.0 , LOCEAN-IPSL (2006)55 !! NEMO/OPA 3,2 , LOCEAN-IPSL (2009) 55 56 !! $Id$ 56 57 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 174 175 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & 175 176 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 176 177 ! 177 178 END SUBROUTINE dyn_vor 178 179 … … 206 207 INTEGER , INTENT(in ) :: kt ! ocean time-step index 207 208 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; 208 !! =nrvm (relative vorticity or metric)209 ! ! =nrvm (relative vorticity or metric) 209 210 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend 210 211 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend … … 222 223 ENDIF 223 224 224 ! Local constant initialization 225 zfact2 = 0.5 * 0.5 225 zfact2 = 0.5 * 0.5 ! Local constant initialization 226 226 227 227 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) … … 229 229 DO jk = 1, jpkm1 ! Horizontal slab 230 230 ! ! =============== 231 ! 231 232 ! Potential vorticity and horizontal fluxes 232 233 ! ----------------------------------------- … … 315 316 INTEGER, INTENT(in) :: kt ! ocean timestep index 316 317 !! 317 INTEGER :: ji, jj, jk ! dummy loop indices318 INTEGER :: ji, jj, jk ! dummy loop indices 318 319 REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! temporary scalars 319 320 REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! " " … … 327 328 ENDIF 328 329 329 ! Local constant initialization 330 zfact1 = 0.5 * 0.25 330 zfact1 = 0.5 * 0.25 ! Local constant initialization 331 331 zfact2 = 0.5 * 0.5 332 332 … … 335 335 DO jk = 1, jpkm1 ! Horizontal slab 336 336 ! ! =============== 337 337 ! 338 338 ! Relative and planetary potential vorticity and horizontal fluxes 339 339 ! ---------------------------------------------------------------- … … 438 438 ENDIF 439 439 440 ! Local constant initialization 441 zfact1 = 0.5 * 0.25 440 zfact1 = 0.5 * 0.25 ! Local constant initialization 442 441 443 442 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) … … 445 444 DO jk = 1, jpkm1 ! Horizontal slab 446 445 ! ! =============== 446 ! 447 447 ! Potential vorticity and horizontal fluxes 448 448 ! ----------------------------------------- … … 465 465 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 466 466 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 467 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &467 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 468 468 & ) 469 469 END DO 470 470 END DO 471 471 END SELECT 472 472 ! 473 473 IF( ln_sco ) THEN 474 474 DO jj = 1, jpj ! caution: don't use (:,:) for this loop … … 487 487 END DO 488 488 ENDIF 489 489 ! 490 490 ! Compute and add the vorticity term trend 491 491 ! ---------------------------------------- … … 514 514 !! 515 515 !! ** Method : Trend evaluated using now fields (centered in time) 516 !! and the Arakawa and Lamb (19 XX) flux form formulation : conserves516 !! and the Arakawa and Lamb (1980) flux form formulation : conserves 517 517 !! both the horizontal kinetic energy and the potential enstrophy 518 !! when horizontal divergence is zero. 519 !! The trend of the vorticity term is given by: 520 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 521 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 522 !! Add this trend to the general momentum trend (ua,va): 523 !! (ua,va) = (ua,va) + ( voru , vorv ) 518 !! when horizontal divergence is zero (see the NEMO documentation) 519 !! Add this trend to the general momentum trend (ua,va). 524 520 !! 525 521 !! ** Action : - Update (ua,va) with the now vorticity term trend … … 531 527 INTEGER , INTENT(in ) :: kt ! ocean time-step index 532 528 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; 533 !! =nrvm (relative vorticity or metric)529 ! ! =nrvm (relative vorticity or metric) 534 530 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend 535 531 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 536 532 !! 537 INTEGER :: ji, jj, jk! dummy loop indices533 INTEGER :: ji, jj, jk ! dummy loop indices 538 534 REAL(wp) :: zfac12, zua, zva ! temporary scalars 539 535 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! temporary 2D workspace 540 536 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse ! temporary 3D workspace 537 #if defined key_vvl 538 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3f 539 #else 541 540 REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE :: ze3f 541 #endif 542 542 !!---------------------------------------------------------------------- 543 543 … … 546 546 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 547 547 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 548 548 ENDIF 549 550 IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t) 549 551 DO jk = 1, jpk 550 552 DO jj = 1, jpjm1 … … 559 561 ENDIF 560 562 561 ! Local constant initialization 562 zfac12 = 1.e0 / 12.e0 563 zfac12 = 1.e0 / 12.e0 ! Local constant initialization 563 564 564 565 … … 571 572 ! ----------------------------------------- 572 573 SELECT CASE( kvor ) ! vorticity considered 573 CASE ( 1 ) ; zwz(:,:) = ff(:,:) * ze3f(:,:,jk) ! planetary vorticity (Coriolis) 574 CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) ! relative vorticity 574 CASE ( 1 ) ! planetary vorticity (Coriolis) 575 zwz(:,:) = ff(:,:) * ze3f(:,:,jk) 576 CASE ( 2 ) ! relative vorticity 577 zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 575 578 CASE ( 3 ) ! metric term 576 579 DO jj = 1, jpjm1 … … 581 584 END DO 582 585 END DO 583 CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) ! total (relative + planetary vorticity) 586 CASE ( 4 ) ! total (relative + planetary vorticity) 587 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 584 588 CASE ( 5 ) ! total (coriolis + metric) 585 589 DO jj = 1, jpjm1 … … 588 592 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 589 593 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 590 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &594 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 591 595 & ) * ze3f(ji,jj,jk) 592 596 END DO … … 599 603 ! Compute and add the vorticity term trend 600 604 ! ---------------------------------------- 601 jj =2602 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ;ztsw(1,:) = 0605 jj = 2 606 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 603 607 DO ji = 2, jpi 604 608 ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) … … 636 640 !! 637 641 !! ** Purpose : Control the consistency between cpp options for 638 !! tracer advection schemes642 !! tracer advection schemes 639 643 !!---------------------------------------------------------------------- 640 644 INTEGER :: ioptio ! temporary integer -
trunk/NEMO/OPA_SRC/DYN/wzvmod.F90
r1241 r1438 1 1 MODULE wzvmod 2 !! MODULE sshwzv 2 3 !!============================================================================== 3 !! *** MODULE wzvmod***4 !! Ocean d iagnostic variable :vertical velocity4 !! *** MODULE sshwzv *** 5 !! Ocean dynamics : sea surface height and vertical velocity 5 6 !!============================================================================== 6 !! History : 5.0 ! 90-10 (C. Levy, G. Madec) Original code7 !! 7.0 ! 96-01 (G. Madec) Statement function for e38 !! 8.5 ! 02-07 (G. Madec) Free form, F90 9 !! " ! 07-07 (D. Storkey) Zero zhdiv at open boundary (BDY)10 !! ----------------------------------------------------------------------11 !! wzv : Compute the vertical velocity12 !! ----------------------------------------------------------------------13 !! * Modules used7 !! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code 8 !!---------------------------------------------------------------------- 9 10 !!---------------------------------------------------------------------- 11 !! ssh_wzv : after ssh & now vertical velocity 12 !! ssh_nxt : filter ans swap the ssh arrays 13 !! ssh_rst : read/write ssh restart fields in the ocean restart file 14 !!---------------------------------------------------------------------- 14 15 USE oce ! ocean dynamics and tracers variables 15 16 USE dom_oce ! ocean space and time domain variables 16 17 USE sbc_oce ! surface boundary condition: ocean 17 18 USE domvvl ! Variable volume 19 USE iom ! I/O library 20 USE restart ! only for lrst_oce 18 21 USE in_out_manager ! I/O manager 19 22 USE prtctl ! Print control 20 23 USE phycst 21 USE bdy_oce ! unstructured open boundaries22 24 USE lbclnk ! ocean lateral boundary condition (or mpp link) 23 25 USE obc_par ! open boundary cond. parameter … … 27 29 PRIVATE 28 30 29 !! * Routine accessibility30 PUBLIC wzv ! routine called by step.F90 and inidtr.F9031 PUBLIC ssh_wzv ! called by step.F90 32 PUBLIC ssh_nxt ! called by step.F90 31 33 32 34 !! * Substitutions 33 35 # include "domzgr_substitute.h90" 34 !!---------------------------------------------------------------------- 35 !! OPA 9.0 , LOCEAN-IPSL (2005) 36 # include "vectopt_loop_substitute.h90" 37 38 !!---------------------------------------------------------------------- 39 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 36 40 !! $Id$ 37 41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 40 44 CONTAINS 41 45 42 SUBROUTINE wzv( kt )43 !!---------------------------------------------------------------------- 44 !! *** ROUTINEwzv ***45 !! 46 !! ** Purpose : Compute the now vertical velocity after the array swap47 !! 48 !! ** Method : Using the incompressibility hypothesis, the vertical49 !! velocity is computed by integrating the horizontal divergence50 !! from the bottom to the surface.51 !! The boundary conditions are w=0 at the bottom (no flux) and,52 !! in regid-lid case, w=0 at the sea surface.53 !! 54 !! ** action : wn array : the now vertical velocity55 !! ----------------------------------------------------------------------56 !! * Arguments57 INTEGER, INTENT( in ) :: kt ! ocean time-step index58 59 !! * Local declarations60 INTEGER :: jk ! dummy loop indices61 !! Variable volume62 INTEGER :: ji, jj ! dummy loop indices63 REAL(wp) :: z2dt, zraur ! temporary scalar64 REAL(wp), DIMENSION (jpi,jpj) :: zssha, zun, zvn, zhdiv65 #if defined key_bdy 66 INTEGER :: jgrd, jb! temporary scalars67 #endif 46 SUBROUTINE ssh_wzv( kt ) 47 !!---------------------------------------------------------------------- 48 !! *** ROUTINE ssh_wzv *** 49 !! 50 !! ** Purpose : compute the after ssh (ssha), the now vertical velocity 51 !! and update the now vertical coordinate (lk_vvl=T). 52 !! 53 !! ** Method : - 54 !! 55 !! - Using the incompressibility hypothesis, the vertical 56 !! velocity is computed by integrating the horizontal divergence 57 !! from the bottom to the surface minus the scale factor evolution. 58 !! The boundary conditions are w=0 at the bottom (no flux) and. 59 !! 60 !! ** action : ssha : after sea surface height 61 !! wn : now vertical velocity 62 !! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height 63 !! at u-, v-, f-point s 64 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 65 !!---------------------------------------------------------------------- 66 INTEGER, INTENT(in) :: kt ! time step 67 !! 68 INTEGER :: ji, jj, jk ! dummy loop indices 69 REAL(wp) :: zcoefu, zcoefv, zcoeff ! temporary scalars 70 REAL(wp) :: z2dt, zraur ! temporary scalars 71 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace 68 72 !!---------------------------------------------------------------------- 69 73 70 74 IF( kt == nit000 ) THEN 71 75 IF(lwp) WRITE(numout,*) 72 IF(lwp) WRITE(numout,*) 'wzv : vertical velocity from continuity eq.' 73 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 74 75 ! bottom boundary condition: w=0 (set once for all) 76 wn(:,:,jpk) = 0.e0 77 ENDIF 78 79 IF( lk_vvl ) THEN ! Variable volume 80 ! 81 z2dt = 2. * rdt ! time step: leap-frog 82 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest 83 zraur = 1. / rauw 84 85 ! Vertically integrated quantities 86 ! -------------------------------- 87 zun(:,:) = 0.e0 88 zvn(:,:) = 0.e0 89 ! 90 DO jk = 1, jpkm1 ! Vertically integrated transports (now) 91 zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk) 92 zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 93 END DO 94 95 ! Horizontal divergence of barotropic transports 96 !-------------------------------------------------- 97 zhdiv(:,:) = 0.e0 98 DO jj = 2, jpjm1 99 DO ji = 2, jpim1 ! vector opt. 100 zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun(ji ,jj ) & 101 & - e2u(ji-1,jj ) * zun(ji-1,jj ) & 102 & + e1v(ji ,jj ) * zvn(ji ,jj ) & 103 & - e1v(ji ,jj-1) * zvn(ji ,jj-1) ) & 104 & / ( e1t(ji,jj) * e2t(ji,jj) ) 76 IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' 77 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 78 ! 79 CALL ssh_rst( nit000, 'READ' ) ! read or initialize sshb and sshn 80 ! 81 wn(:,:,jpk) = 0.e0 ! bottom boundary condition: w=0 (set once for all) 82 ! 83 IF( lk_vvl ) THEN ! before and now Sea SSH at u-, v-, f-points (vvl case only) 84 DO jj = 1, jpjm1 85 DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible 86 zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 87 zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 88 zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 89 sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 90 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 91 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 92 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 93 sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) & 94 & + sshb(ji+1,jj) + sshb(ji+1,jj+1) ) 95 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & 96 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 97 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & 98 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 99 sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) & 100 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) 101 END DO 102 END DO 103 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) 104 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) 105 CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) 106 ENDIF 107 ! 108 ENDIF 109 110 ! set time step size (Euler/Leapfrog) 111 z2dt = 2. * rdt 112 IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt 113 114 zraur = 1. / rauw 115 116 ! !------------------------------! 117 ! ! After Sea Surface Height ! 118 ! !------------------------------! 119 zhdiv(:,:) = 0.e0 120 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports 121 zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 122 END DO 123 124 ! ! Sea surface elevation time stepping 125 ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 126 127 #if defined key_obc 128 # if defined key_agrif 129 IF ( Agrif_Root() ) THEN 130 # endif 131 ssha(:,:) = ssha(:,:) * obctmsk(:,:) 132 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) 133 # if defined key_agrif 134 ENDIF 135 # endif 136 #endif 137 138 ! ! Sea Surface Height at u-,v- and f-points (vvl case only) 139 IF( lk_vvl ) THEN ! (required only in key_vvl case) 140 DO jj = 1, jpjm1 141 DO ji = 1, fs_jpim1 ! Vector Opt. 142 sshu_a(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 143 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) & 144 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 145 sshv_a(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 146 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 147 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 148 sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1) & ! Caution : fmask not used 149 & * ( ssha(ji ,jj) + ssha(ji ,jj+1) & 150 & + ssha(ji+1,jj) + ssha(ji+1,jj+1) ) 105 151 END DO 106 152 END DO 107 108 #if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts ) 109 ! open boundaries (div must be zero behind the open boundary) 110 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 111 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east 112 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west 113 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north 114 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south 115 #endif 116 117 #if defined key_bdy 118 jgrd=1 !: tracer grid. 119 DO jb = 1, nblenrim(jgrd) 120 ji = nbi(jb,jgrd) 121 jj = nbj(jb,jgrd) 122 zhdiv(ji,jj) = 0.e0 153 CALL lbc_lnk( sshu_a, 'U', 1. ) ! Boundaries conditions 154 CALL lbc_lnk( sshv_a, 'V', 1. ) 155 CALL lbc_lnk( sshf_a, 'F', 1. ) 156 ENDIF 157 158 ! !------------------------------! 159 ! ! Now Vertical Velocity ! 160 ! !------------------------------! 161 ! ! integrate from the bottom the hor. divergence 162 DO jk = jpkm1, 1, -1 163 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 164 & - ( fse3t_a(:,:,jk) & 165 & - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 166 END DO 167 168 ! !------------------------------! 169 ! ! Update Now Vertical coord. ! 170 ! !------------------------------! 171 IF( lk_vvl ) THEN ! only in vvl case) 172 ! ! now local depth and scale factors (stored in fse3. arrays) 173 DO jk = 1, jpkm1 174 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! depths 175 fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 176 fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 177 ! 178 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors 179 fse3u (:,:,jk) = fse3u_n (:,:,jk) 180 fse3v (:,:,jk) = fse3v_n (:,:,jk) 181 fse3f (:,:,jk) = fse3f_n (:,:,jk) 182 fse3w (:,:,jk) = fse3w_n (:,:,jk) 183 fse3uw(:,:,jk) = fse3uw_n(:,:,jk) 184 fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 123 185 END DO 124 #endif 125 126 CALL lbc_lnk( zhdiv, 'T', 1. ) 127 128 ! Sea surface elevation time stepping 129 ! ----------------------------------- 130 zssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1) 131 132 ! Vertical velocity computed from bottom 133 ! -------------------------------------- 134 DO jk = jpkm1, 1, -1 135 wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) & 136 & - ( zssha(:,:) - sshb(:,:) ) * fsve3t(:,:,jk) * mut(:,:,jk) / z2dt 137 END DO 138 139 ELSE ! Fixed volume 140 141 ! Vertical velocity computed from bottom 142 ! -------------------------------------- 143 DO jk = jpkm1, 1, -1 144 wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) 145 END DO 146 147 ENDIF 148 149 IF(ln_ctl) CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 - : ', mask1=wn) 150 151 END SUBROUTINE wzv 186 ! ! ocean depth (at u- and v-points) 187 hu(:,:) = hu_0(:,:) + sshu_n(:,:) 188 hv(:,:) = hv_0(:,:) + sshv_n(:,:) 189 ! ! masked inverse of the ocean depth (at u- and v-points) 190 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) ) 191 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) ) 192 193 ENDIF 194 ! 195 END SUBROUTINE ssh_wzv 196 197 198 SUBROUTINE ssh_nxt( kt ) 199 !!---------------------------------------------------------------------- 200 !! *** ROUTINE ssh_nxt *** 201 !! 202 !! ** Purpose : achieve the sea surface height time stepping by 203 !! applying Asselin time filter and swapping the arrays 204 !! ssha already computed in ssh_wzv 205 !! 206 !! ** Method : - apply Asselin time fiter to now ssh and swap : 207 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 208 !! sshn = ssha 209 !! 210 !! ** action : - sshb, sshn : before & now sea surface height 211 !! ready for the next time step 212 !!---------------------------------------------------------------------- 213 INTEGER, INTENT( in ) :: kt ! ocean time-step index 214 !! 215 INTEGER :: ji, jj ! dummy loop indices 216 !!---------------------------------------------------------------------- 217 218 IF( kt == nit000 ) THEN 219 IF(lwp) WRITE(numout,*) 220 IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)' 221 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 222 ENDIF 223 224 ! Time filter and swap of the ssh 225 ! ------------------------------- 226 227 IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points 228 ! ! ---------------------- ! 229 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 230 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) 231 sshu_n(:,:) = sshu_a(:,:) 232 sshv_n(:,:) = sshv_a(:,:) 233 sshf_n(:,:) = sshf_a(:,:) 234 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 235 DO jj = 1, jpj 236 DO ji = 1, jpi ! before <-- now filtered 237 sshb (ji,jj) = sshn(ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) ) 238 sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 239 sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) 240 sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) ) 241 sshn (ji,jj) = ssha (ji,jj) ! now <-- after 242 sshu_n(ji,jj) = sshu_a(ji,jj) 243 sshv_n(ji,jj) = sshv_a(ji,jj) 244 sshf_n(ji,jj) = sshf_a(ji,jj) 245 END DO 246 END DO 247 ENDIF 248 ! 249 ELSE ! fixed levels : ssh at t-point only 250 ! ! ------------ ! 251 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 252 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 253 ! 254 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 255 DO jj = 1, jpj 256 DO ji = 1, jpi ! before <-- now filtered 257 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 258 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 259 END DO 260 END DO 261 ENDIF 262 ! 263 ENDIF 264 265 ! ! write filtered free surface arrays in restart file 266 IF( lrst_oce ) CALL ssh_rst( kt, 'WRITE' ) 267 ! 268 IF(ln_ctl) CALL prt_ctl(tab2d_1=sshb , clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 269 ! 270 END SUBROUTINE ssh_nxt 271 272 273 SUBROUTINE ssh_rst( kt, cdrw ) 274 !!--------------------------------------------------------------------- 275 !! *** ROUTINE ssh_rst *** 276 !! 277 !! ** Purpose : Read or write Sea Surface Height arrays in restart file 278 !! 279 !! ** action : - cdrw = READ : sshb, sshn read in ocean restart file 280 !! - cdrw = WRITE : sshb, sshn written in ocean restart file 281 !!---------------------------------------------------------------------- 282 INTEGER , INTENT(in) :: kt ! ocean time-step 283 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 284 !!---------------------------------------------------------------------- 285 ! 286 IF( TRIM(cdrw) == 'READ' ) THEN 287 IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 288 CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb(:,:) ) 289 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn(:,:) ) 290 IF( neuler == 0 ) sshb(:,:) = sshn(:,:) 291 ELSE 292 IF( nn_rstssh == 1 ) THEN 293 sshb(:,:) = 0.e0 294 sshn(:,:) = 0.e0 295 ENDIF 296 ENDIF 297 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 298 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb(:,:) ) 299 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn(:,:) ) 300 ENDIF 301 ! 302 END SUBROUTINE ssh_rst 152 303 153 304 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
