Changeset 1361 for branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynnxt.F90
- Timestamp:
- 2009-03-31T18:26:41+02:00 (15 years ago)
- File:
-
- 1 edited
-
branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynnxt.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynnxt.F90
r1152 r1361 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.1 ! 2009-02 (G. Madec) re-introduce the vvl option 18 !!---------------------------------------------------------------------- 19 7 20 !!---------------------------------------------------------------------- 8 21 !! dyn_nxt : update the horizontal velocity from the momentum trend … … 34 47 # include "domzgr_substitute.h90" 35 48 !!---------------------------------------------------------------------- 49 !! OPA 9.0 , LOCEAN-IPSL (2005) 50 !! $Id$ 51 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 52 !!---------------------------------------------------------------------- 36 53 37 54 CONTAINS … … 59 76 !! - Update un,vn arrays, the now horizontal velocity 60 77 !! 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 78 !!---------------------------------------------------------------------- 73 !! * Arguments74 79 INTEGER, INTENT( in ) :: kt ! ocean time-step index 75 76 !! * Local declarations 77 INTEGER :: ji, jj, jk ! dummy loop indices 80 !! 81 INTEGER :: jk ! dummy loop indices 78 82 REAL(wp) :: z2dt ! temporary scalar 79 REAL(wp) :: zsshun1, zsshvn1 ! temporary scalar80 !! Variable volume81 REAL(wp), DIMENSION(jpi,jpj) :: & ! 2D workspace82 zsshub, zsshun, zsshua, &83 zsshvb, zsshvn, zsshva84 REAL(wp), DIMENSION(jpi,jpj,jpk) :: &85 zfse3ub, zfse3un, zfse3ua, & ! 3D workspace86 zfse3vb, zfse3vn, zfse3va87 !!----------------------------------------------------------------------88 !! OPA 9.0 , LOCEAN-IPSL (2005)89 !! $Id$90 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt91 83 !!---------------------------------------------------------------------- 92 84 … … 102 94 103 95 !! Explicit physics with thickness weighted updates 104 IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN105 106 ! Sea surface elevation time stepping107 ! -----------------------------------108 !109 DO jj = 1, jpjm1110 DO ji = 1,jpim1111 112 ! Sea Surface Height at u-point before113 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 before118 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 after123 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 after128 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 DO133 END DO134 135 ! Boundaries conditions136 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 step140 ! -------------------------------------------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 step145 ! ----------------------------------------------146 IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN147 CALL dom_vvl_sf_ini( 'U', zfse3un ) ; CALL dom_vvl_sf_ini( 'V', zfse3vn )148 ELSE149 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 ENDIF153 154 ! Thickness weighting155 ! -------------------156 DO jk = 1, jpkm1157 DO jj = 1, jpj158 DO ji = 1, jpi159 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 DO168 END DO169 END DO170 171 ENDIF172 96 173 97 ! Lateral boundary conditions on ( ua, va ) … … 175 99 CALL lbc_lnk( va, 'V', -1. ) 176 100 177 ! ! =============== 178 DO jk = 1, jpkm1 ! Horizontal slab 179 ! ! =============== 180 ! Next velocity 181 ! ------------- 101 ! Next velocity 102 ! ------------- 182 103 #if defined key_dynspg_flt 183 ! Leap-frog time stepping already done in dynspg.F routine104 ! Leap-frog time stepping already done in dynspg_flt.F routine 184 105 #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 106 IF( lk_vvl ) THEN ! Varying levels 107 !RB_vvl scale factors are wrong at this point 108 DO jk = 1, jpkm1 109 ua(ji,jj,jk) = ( ub(:,:,jk) * fse3u(:,:,jk) & 110 & + z2dt * ua(:,:,jk) * fse3u(:,:,jk) ) & 111 & / fse3u(:,:,jk) * umask(:,:,jk) 112 va(ji,jj,jk) = ( vb(:,:,jk) * fse3v(:,:,jk) & 113 & + z2dt * va(:,:,jk) * fse3v(:,:,jk) ) & 114 & / fse3v(:,:,jk) * vmask(:,:,jk) 115 END DO 116 ELSE 117 DO jk = 1, jpkm1 118 ! Leap-frog time stepping 119 ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 120 va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 121 END DO 122 ENDIF 210 123 211 124 # if defined key_obc 212 ! ! ===============213 END DO ! End of slab214 ! ! ===============215 125 ! Update (ua,va) along open boundaries (only in the rigid-lid case) 216 126 CALL obc_dyn( kt ) … … 235 145 ENDIF 236 146 237 ! ! ===============238 DO jk = 1, jpkm1 ! Horizontal slab239 ! ! ===============240 147 # elif defined key_bdy 241 ! ! ===============242 END DO ! End of slab243 ! ! ===============244 148 ! Update (ua,va) along open boundaries (for exp or ts options). 245 149 IF ( lk_dynspg_exp .or. lk_dynspg_ts ) THEN … … 279 183 ENDIF 280 184 281 ! ! ===============282 DO jk = 1, jpkm1 ! Horizontal slab283 ! ! ===============284 185 # endif 285 186 # 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 187 CALL Agrif_dyn( kt ) 291 ! ! ===============292 DO jk = 1, jpkm1 ! Horizontal slab293 ! ! ===============294 188 # endif 295 189 #endif 296 190 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 303 DO jj = 1, jpj 304 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) 313 END DO 314 END DO 315 ELSE ! Fixed levels 316 DO jj = 1, jpj 317 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) 321 un(ji,jj,jk) = ua(ji,jj,jk) 322 vn(ji,jj,jk) = va(ji,jj,jk) 323 END DO 324 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 191 ! Time filter and swap of dynamics arrays 192 ! ------------------------------------------ 193 IF( neuler == 0 .AND. kt == nit000 ) THEN 194 DO jk = 1, jpkm1 195 ub(:,:,jk) = un(:,:,jk) 196 vb(:,:,jk) = vn(:,:,jk) 197 un(:,:,jk) = ua(:,:,jk) 198 vn(:,:,jk) = va(:,:,jk) 199 END DO 200 ELSE 201 IF( lk_vvl ) THEN ! Varying levels 202 ! Not done 203 ELSE ! Fixed levels 204 !RB_vvl : should be done as in tranxt ? 205 DO jk = 1, jpkm1 ! filter applied on velocities 206 ub(:,:,jk) = atfp * ( ub(:,:,jk) + ua(:,:,jk) ) + atfp1 * un(:,:,jk) 207 vb(:,:,jk) = atfp * ( vb(:,:,jk) + va(:,:,jk) ) + atfp1 * vn(:,:,jk) 208 un(:,:,jk) = ua(:,:,jk) 209 vn(:,:,jk) = va(:,:,jk) 210 END DO 355 211 ENDIF 356 ! ! =============== 357 END DO ! End of slab 358 ! ! =============== 212 ENDIF 359 213 360 214 IF(ln_ctl) THEN
Note: See TracChangeset
for help on using the changeset viewer.
