Changeset 1601 for trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
- Timestamp:
- 2009-08-11T12:09:19+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r1556 r1601 18 18 !! 'key_dynspg_flt' filtered free surface 19 19 !!---------------------------------------------------------------------- 20 !! dyn_spg_flt : update the momentum trend with the surface pressure 21 !! gradient in the filtered free surface case with 22 !! vector optimization 20 !! dyn_spg_flt : update the momentum trend with the surface pressure gradient in the filtered free surface case 23 21 !! flt_rst : read/write the time-splitting restart fields in the ocean restart file 24 22 !!---------------------------------------------------------------------- … … 31 29 USE phycst ! physical constants 32 30 USE domvvl ! variable volume 31 USE solmat ! matrix construction for elliptic solvers 33 32 USE solver ! solver initialization 34 33 USE solpcg ! preconditionned conjugate gradient solver … … 44 43 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 45 44 USE prtctl ! Print control 46 USE solmat ! matrix construction for elliptic solvers47 45 USE agrif_opa_interp 48 46 USE iom … … 76 74 !! ** Method : Filtered free surface formulation. The surface 77 75 !! pressure gradient is given by: 78 !! spgu = 1/rau0 d/dx(ps) = 1/e1u di( sshn + rnubtda )79 !! spgv = 1/rau0 d/dy(ps) = 1/e2v dj( sshn + rnubtda )76 !! spgu = 1/rau0 d/dx(ps) = 1/e1u di( sshn + btda ) 77 !! spgv = 1/rau0 d/dy(ps) = 1/e2v dj( sshn + btda ) 80 78 !! where sshn is the free surface elevation and btda is the after 81 79 !! time derivative of the free surface elevation … … 106 104 USE oce, ONLY : zvb => sa ! ta used as workspace 107 105 !! 108 INTEGER, INTENT( in ):: kt ! ocean time-step index109 INTEGER, INTENT( out) :: kindic ! solver convergence flag (<0 if not converge)106 INTEGER, INTENT(in ) :: kt ! ocean time-step index 107 INTEGER, INTENT( out) :: kindic ! solver convergence flag (<0 if not converge) 110 108 !! 111 INTEGER :: ji, jj, jk 112 REAL(wp) :: z2dt, z2dtg, zraur , znugdt! temporary scalars113 REAL(wp) :: z nurau, zgcb, zbtd ! " "114 REAL(wp) :: ztdgu, ztdgv ! " "109 INTEGER :: ji, jj, jk ! dummy loop indices 110 REAL(wp) :: z2dt, z2dtg, zraur ! temporary scalars 111 REAL(wp) :: zgcb, zbtd ! - - 112 REAL(wp) :: ztdgu, ztdgv ! - - 115 113 !!---------------------------------------------------------------------- 116 114 ! … … 127 125 ! read filtered free surface arrays in restart file 128 126 ! when using agrif, sshn, gcx have to be read in istate 129 IF (.NOT. lk_agrif) CALL flt_rst( nit000, 'READ' )! read or initialize the following fields:127 IF(.NOT. lk_agrif) CALL flt_rst( nit000, 'READ' ) ! read or initialize the following fields: 130 128 ! ! gcx, gcxb 131 129 ENDIF 132 130 133 131 ! Local constant initialization 134 z2dt = 2. * rdt ! time step: leap-frog135 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt! time step: Euler if restart from rest136 IF( neuler == 0 .AND. kt == nit000+1 ) CALL sol_mat(kt)132 z2dt = 2. * rdt ! time step: leap-frog 133 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest 134 IF( neuler == 0 .AND. kt == nit000+1 ) CALL sol_mat( kt ) 137 135 z2dtg = grav * z2dt 138 136 zraur = 1. / rauw 139 znugdt = rnu * grav * z2dt140 znurau = znugdt * zraur141 137 142 138 !! Explicit physics with thickness weighted updates … … 237 233 END DO 238 234 END DO 239 240 ! Boundary conditions on (spgu,spgv) 241 CALL lbc_lnk( spgu, 'U', -1. ) 235 CALL lbc_lnk( spgu, 'U', -1. ) ! lateral boundary conditions 242 236 CALL lbc_lnk( spgv, 'V', -1. ) 243 237 … … 245 239 246 240 ! Right hand side of the elliptic equation and first guess 247 ! -------------------------------------------------------- ---241 ! -------------------------------------------------------- 248 242 DO jj = 2, jpjm1 249 243 DO ji = fs_2, fs_jpim1 ! vector opt. … … 259 253 END DO 260 254 ! applied the lateral boundary conditions 261 IF( n solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 )CALL lbc_lnk_e( gcb, c_solver_pt, 1. )255 IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 ) CALL lbc_lnk_e( gcb, c_solver_pt, 1. ) 262 256 263 257 #if defined key_agrif … … 265 259 ! add contribution of gradient of after barotropic transport divergence 266 260 IF( nbondi == -1 .OR. nbondi == 2 ) gcb(3 ,:) = & 267 & gcb(3 ,:) - z nugdt* z2dt * laplacu(2 ,:) * gcdprc(3 ,:) * hu(2 ,:) * e2u(2 ,:)261 & gcb(3 ,:) - z2dtg * z2dt * laplacu(2 ,:) * gcdprc(3 ,:) * hu(2 ,:) * e2u(2 ,:) 268 262 IF( nbondi == 1 .OR. nbondi == 2 ) gcb(nlci-2,:) = & 269 & gcb(nlci-2,:) + z nugdt* z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)263 & gcb(nlci-2,:) + z2dtg * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:) 270 264 IF( nbondj == -1 .OR. nbondj == 2 ) gcb(: ,3) = & 271 & gcb(:,3 ) - z nugdt* z2dt * laplacv(:,2 ) * gcdprc(:,3 ) * hv(:,2 ) * e1v(:,2 )265 & gcb(:,3 ) - z2dtg * z2dt * laplacv(:,2 ) * gcdprc(:,3 ) * hv(:,2 ) * e1v(:,2 ) 272 266 IF( nbondj == 1 .OR. nbondj == 2 ) gcb(:,nlcj-2) = & 273 & gcb(:,nlcj-2) + z nugdt* z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)267 & gcb(:,nlcj-2) + z2dtg * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2) 274 268 ENDIF 275 269 #endif … … 298 292 kindic = 0 299 293 IF( ncut == 0 ) THEN 300 IF( nsolv == 1 ) THEN ! diagonal preconditioned conjuguate gradient 301 CALL sol_pcg( kindic ) 302 ELSEIF( nsolv == 2 ) THEN ! successive-over-relaxation 303 CALL sol_sor( kindic ) 304 ELSE ! e r r o r in nsolv namelist parameter 305 WRITE(ctmp1,*) ' ~~~~~~~~~~~ not = ', nsolv 306 CALL ctl_stop( ' dyn_spg_flt : e r r o r, nsolv = 1 or 2', ctmp1 ) 294 IF ( nn_solv == 1 ) THEN ; CALL sol_pcg( kindic ) ! diagonal preconditioned conjuguate gradient 295 ELSEIF( nn_solv == 2 ) THEN ; CALL sol_sor( kindic ) ! successive-over-relaxation 307 296 ENDIF 308 297 ENDIF … … 313 302 DO ji = fs_2, fs_jpim1 ! vector opt. 314 303 ! trend of Transport divergence gradient 315 ztdgu = z nugdt* (gcx(ji+1,jj ) - gcx(ji,jj) ) / e1u(ji,jj)316 ztdgv = z nugdt* (gcx(ji ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)304 ztdgu = z2dtg * (gcx(ji+1,jj ) - gcx(ji,jj) ) / e1u(ji,jj) 305 ztdgv = z2dtg * (gcx(ji ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj) 317 306 ! multiplied by z2dt 318 307 #if defined key_obc … … 336 325 IF( .NOT. Agrif_Root() ) THEN 337 326 ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface 338 IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2 ,:) = z nugdt* z2dt * laplacu(2 ,:) * umask(2 ,:,1)339 IF( nbondi == 1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z nugdt* z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)340 IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2 ) = z nugdt* z2dt * laplacv(:,2 ) * vmask(: ,2,1)341 IF( nbondj == 1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z nugdt* z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)327 IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2 ,:) = z2dtg * z2dt * laplacu(2 ,:) * umask(2 ,:,1) 328 IF( nbondi == 1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z2dtg * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1) 329 IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2 ) = z2dtg * z2dt * laplacv(:,2 ) * vmask(: ,2,1) 330 IF( nbondj == 1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z2dtg * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1) 342 331 ENDIF 343 332 #endif
Note: See TracChangeset
for help on using the changeset viewer.
