Changeset 1438 for trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
- Timestamp:
- 2009-05-11T16:34:47+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
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 !
Note: See TracChangeset
for help on using the changeset viewer.
