Changeset 1390
- Timestamp:
- 2009-04-08T16:44:58+02:00 (15 years ago)
- Location:
- branches/dev_004_VVL/NEMO/OPA_SRC/DYN
- Files:
-
- 4 edited
-
dynnxt.F90 (modified) (1 diff)
-
dynspg_exp.F90 (modified) (6 diffs)
-
dynspg_flt.F90 (modified) (7 diffs)
-
wzvmod.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynnxt.F90
r1382 r1390 110 110 IF( lk_vvl ) THEN ! Varying levels 111 111 DO jk = 1, jpkm1 112 ua( ji,jj,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) &112 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 113 113 & + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 114 114 & / fse3u_a(:,:,jk) * umask(:,:,jk) 115 va( ji,jj,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) &115 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 116 116 & + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 117 117 & / fse3v_a(:,:,jk) * vmask(:,:,jk) -
branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r1146 r1390 11 11 !! pressure gradient in the free surface constant 12 12 !! volume case with vector optimization 13 !! exp_rst : read/write the explicit restart fields in the ocean restart file14 13 !!---------------------------------------------------------------------- 15 14 !! * Modules used … … 33 32 !! * Accessibility 34 33 PUBLIC dyn_spg_exp ! routine called by step.F90 35 PUBLIC exp_rst ! routine called by istate.F9036 34 37 35 !! * Substitutions … … 62 60 !! -1- Compute the now surface pressure gradient 63 61 !! -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 62 !! 80 63 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend … … 91 74 !! * Local declarations 92 75 INTEGER :: ji, jj, jk ! dummy loop indices 93 REAL(wp) :: z2dt, zraur, zssha ! temporary scalars94 REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary arrays95 & zhdiv96 76 !!---------------------------------------------------------------------- 97 77 … … 104 84 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 105 85 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 86 ENDIF 111 87 112 IF( .NOT. lk_vvl ) THEN 88 ! 0. Initialization 89 ! ----------------- 90 ! read or estimate sea surface height and vertically integrated velocities 91 IF( lk_obc ) CALL obc_dta_bt( kt, 0 ) 113 92 114 ! 0. Initialization 115 ! ----------------- 116 ! read or estimate sea surface height and vertically integrated velocities 117 IF( lk_obc ) CALL obc_dta_bt( kt, 0 ) 118 z2dt = 2. * rdt ! time step: leap-frog 119 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest 120 zraur = 1.e0 / rauw 93 ! 1. Surface pressure gradient (now) 94 ! ---------------------------- 95 DO jj = 2, jpjm1 96 DO ji = fs_2, fs_jpim1 ! vector opt. 97 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj) 98 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 99 END DO 100 END DO 121 101 122 ! 1. Surface pressure gradient (now) 123 ! ---------------------------- 102 ! 2. Add the surface pressure trend to the general trend 103 ! ------------------------------------------------------ 104 DO jk = 1, jpkm1 124 105 DO jj = 2, jpjm1 125 106 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)107 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 108 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 128 109 END DO 129 110 END DO 111 END DO 130 112 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 113 END SUBROUTINE dyn_spg_exp 226 114 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 115 #else 256 116 !!---------------------------------------------------------------------- … … 261 121 WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt 262 122 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 123 #endif 269 124 270 125 !!====================================================================== 271 126 END MODULE dynspg_exp -
branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r1389 r1390 81 81 !! where sshn is the free surface elevation and btda is the after 82 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- 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 !! … … 123 115 & znurau, zgcb, zbtd, & ! " " 124 116 & ztdgu, ztdgv ! " " 125 !! Variable volume126 REAL(wp), DIMENSION(jpi,jpj) :: &127 zsshub, zsshua, zsshvb, zsshva, zssha ! 2D workspace128 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & ! 3D workspace129 zfse3ub, zfse3ua, zfse3vb, zfse3va ! 3D workspace130 117 !!---------------------------------------------------------------------- 131 118 ! … … 182 169 END DO 183 170 END DO 184 185 ! Add the surface pressure trend to the general trend 171 ! 172 ! add the surface pressure trend to the general trend and 173 ! evaluate the masked next velocity (effect of the additional force not included) 186 174 DO jk = 1, jpkm1 187 175 DO jj = 2, jpjm1 188 176 DO ji = fs_2, fs_jpim1 ! vector opt. 189 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)190 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)177 ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) ) ) * umask(ji,jj,jk) 178 va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) ) ) * vmask(ji,jj,jk) 191 179 END DO 192 180 END DO 193 181 END DO 194 195 ! Evaluate the masked next velocity (effect of the additional force not included) 196 DO jk = 1, jpkm1 197 DO jj = 2, jpjm1 198 DO ji = fs_2, fs_jpim1 ! vector opt. 199 ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk) 200 va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk) 201 END DO 202 END DO 203 END DO 204 182 ! 205 183 ENDIF 206 184 … … 266 244 CALL lbc_lnk( spgv, 'V', -1. ) 267 245 268 IF( lk_vvl ) CALL sol_mat( kt ) 246 IF( lk_vvl ) CALL sol_mat( kt ) ! build the matrix at kt (vvl case only) 269 247 270 248 ! Right hand side of the elliptic equation and first guess … … 302 280 ! Relative precision (computation on one processor) 303 281 ! ------------------ 304 rnorme =0. 282 rnorme =0.e0 305 283 rnorme = SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) ) 306 284 IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain … … 368 346 ENDIF 369 347 #endif 370 ! 7.Add the trends multiplied by z2dt to the after velocity371 ! ------------------------------------------------------- ----348 ! Add the trends multiplied by z2dt to the after velocity 349 ! ------------------------------------------------------- 372 350 ! ( c a u t i o n : (ua,va) here are the after velocity not the 373 351 ! trend, the leap-frog time stepping will not -
branches/dev_004_VVL/NEMO/OPA_SRC/DYN/wzvmod.F90
r1389 r1390 154 154 ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 155 155 156 #if defined key_obc 157 # if defined key_agrif 158 IF ( Agrif_Root() ) THEN 159 # endif 160 ssha(:,:) = ssha(:,:)*obctmsk(:,:) 161 CALL lbc_lnk(ssha,'T',1.) ! absolutly compulsory !! (jmm) 162 # if defined key_agrif 163 ENDIF 164 # endif 165 #endif 166 156 167 ! ! Sea Surface Height at u-,v- and f-points (vvl case only) 157 168 IF( lk_vvl ) THEN ! (required only in key_vvl case)
Note: See TracChangeset
for help on using the changeset viewer.
