Changeset 508 for trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
- Timestamp:
- 2006-10-03T17:58:55+02:00 (18 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
r474 r508 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 #if ( defined key_dynspg_flt && ! defined key_mpp_omp ) || defined key_esopa 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 !!---------------------------------------------------------------------- 14 #if ( defined key_dynspg_flt && ! defined key_mpp_omp ) || defined key_esopa 7 15 !!---------------------------------------------------------------------- 8 16 !! 'key_dynspg_flt' filtered free surface 9 17 !! NOT 'key_mpp_omp' k-j-i loop (vector opt.) 18 !!---------------------------------------------------------------------- 10 19 !!---------------------------------------------------------------------- 11 20 !! dyn_spg_flt : update the momentum trend with the surface pressure 12 21 !! gradient in the filtered free surface case with 13 22 !! vector optimization 14 !! ----------------------------------------------------------------------15 !! * Modules used23 !! flt_rst : read/write the time-splitting restart fields in the ocean restart file 24 !!---------------------------------------------------------------------- 16 25 USE oce ! ocean dynamics and tracers 17 26 USE dom_oce ! ocean space and time domain … … 21 30 USE flxrnf ! ocean runoffs 22 31 USE sol_oce ! ocean elliptic solver 32 USE solver ! solver initialization 23 33 USE solpcg ! preconditionned conjugate gradient solver 24 34 USE solsor ! Successive Over-relaxation solver … … 32 42 USE cla_dynspg ! cross land advection 33 43 USE prtctl ! Print control 34 USE in_out_manager ! I/O manager35 44 USE solmat ! matrix construction for elliptic solvers 36 45 USE agrif_opa_interp 46 USE in_out_manager ! I/O manager 47 USE iom 48 USE restart ! only for lrst_oce 37 49 38 50 IMPLICIT NONE 39 51 PRIVATE 40 52 41 !! * Accessibility42 53 PUBLIC dyn_spg_flt ! routine called by step.F90 43 54 … … 48 59 !! OPA 9.0 , LOCEAN-IPSL (2005) 49 60 !! $Header$ 50 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt61 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 51 62 !!---------------------------------------------------------------------- 52 63 … … 96 107 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 97 108 !! 98 !! References : 99 !! Roullet and Madec 1999, JGR. 100 !! 101 !! History : 102 !! ! 98-05 (G. Roullet) Original code 103 !! ! 98-10 (G. Madec, M. Imbard) release 8.2 104 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 105 !! ! 02-11 (C. Talandier, A-M Treguier) Open boundaries 106 !! 9.0 ! 04-08 (C. Talandier) New trends organization 107 !! " ! 05-11 (V. Garnier) Surface pressure gradient organization 109 !! References : Roullet and Madec 1999, JGR. 108 110 !!--------------------------------------------------------------------- 109 !! * Arguments110 111 INTEGER, INTENT( in ) :: kt ! ocean time-step index 111 INTEGER, INTENT( out ) :: kindic ! solver convergence flag 112 ! if the solver doesn t converge 113 ! the flag is < 0 114 !! * Local declarations 115 INTEGER :: ji, jj, jk ! dummy loop indices 116 REAL(wp) :: & 117 z2dt, z2dtg, zraur, znugdt, & ! temporary scalars 118 znurau, zssha, zgcb, zbtd, & ! " " 119 ztdgu, ztdgv ! " " 112 INTEGER, INTENT( out ) :: kindic ! solver convergence flag (<0 if not converge) 113 !! 114 INTEGER :: ji, jj, jk ! dummy loop indices 115 REAL(wp) :: z2dt, z2dtg, zraur, znugdt, & ! temporary scalars 116 & znurau, zssha, zgcb, zbtd, & ! " " 117 & ztdgu, ztdgv ! " " 120 118 !!---------------------------------------------------------------------- 121 119 ! 122 120 IF( kt == nit000 ) THEN 123 121 IF(lwp) WRITE(numout,*) … … 128 126 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 129 127 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 128 CALL solver_init( nit000 ) ! Elliptic solver initialisation 129 130 ! read filtered free surface arrays in restart file 131 CALL flt_rst( nit000, 'READ' ) ! read or initialize the following fields: 132 ! ! gcx, gcxb, sshb, sshn 130 133 ENDIF 131 134 … … 168 171 169 172 #if defined key_obc 170 ! Update velocities on each open boundary with the radiation algorithm 171 CALL obc_dyn( kt ) 172 ! Correction of the barotropic componant velocity to control the volume of the system 173 CALL obc_vol( kt ) 173 CALL obc_dyn( kt ) ! Update velocities on each open boundary with the radiation algorithm 174 CALL obc_vol( kt ) ! Correction of the barotropic componant velocity to control the volume of the system 174 175 #endif 175 176 #if defined key_agrif 176 ! Update velocities on each coarse/fine interfaces 177 178 CALL Agrif_dyn( kt ) 177 CALL Agrif_dyn( kt ) ! Update velocities on each coarse/fine interfaces 179 178 #endif 180 179 #if defined key_orca_r2 … … 243 242 IF( .NOT. AGRIF_ROOT() ) THEN 244 243 ! add contribution of gradient of after barotropic transport divergence 245 IF( (nbondi == -1) .OR. (nbondi == 2) ) gcb(3,:) = gcb(3,:)&246 & -znugdt * z2dt*laplacu(2,:)*gcdprc(3,:)*hu(2,:)*e2u(2,:)247 IF( (nbondi == 1) .OR. (nbondi == 2) ) gcb(nlci-2,:) = gcb(nlci-2,:)&248 & +znugdt * z2dt*laplacu(nlci-2,:)*gcdprc(nlci-2,:)*hu(nlci-2,:)*e2u(nlci-2,:)249 IF( (nbondj == -1) .OR. (nbondj == 2) ) gcb(:,3) = gcb(:,3)&250 & -znugdt * z2dt*laplacv(:,2)*gcdprc(:,3)*hv(:,2)*e1v(:,2)251 IF( (nbondj == 1) .OR. (nbondj == 2) ) gcb(:,nlcj-2) = gcb(:,nlcj-2)&252 & +znugdt * z2dt*laplacv(:,nlcj-2)*gcdprc(:,nlcj-2)*hv(:,nlcj-2)*e1v(:,nlcj-2)244 IF( nbondi == -1 .OR. nbondi == 2 ) gcb(3 ,:) = & 245 & gcb(3 ,:) - znugdt * z2dt * laplacu(2 ,:) * gcdprc(3 ,:) * hu(2 ,:) * e2u(2 ,:) 246 IF( nbondi == 1 .OR. nbondi == 2 ) gcb(nlci-2,:) = & 247 & gcb(nlci-2,:) + znugdt * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:) 248 IF( nbondj == -1 .OR. nbondj == 2 ) gcb(: ,3) = & 249 & gcb(:,3 ) - znugdt * z2dt * laplacv(:,2 ) * gcdprc(:,3 ) * hv(:,2 ) * e1v(:,2 ) 250 IF( nbondj == 1 .OR. nbondj == 2 ) gcb(:,nlcj-2) = & 251 & gcb(:,nlcj-2) + znugdt * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2) 253 252 ENDIF 254 253 #endif … … 263 262 epsr = eps * eps * rnorme 264 263 ncut = 0 265 ! if rnorme is 0, the solution is 0, the solver is n't called264 ! if rnorme is 0, the solution is 0, the solver is not called 266 265 IF( rnorme == 0.e0 ) THEN 267 266 gcx(:,:) = 0.e0 … … 313 312 IF( .NOT. Agrif_Root() ) THEN 314 313 ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface 315 IF( (nbondi == -1) .OR. (nbondi == 2) ) spgu(2,:) = znugdt * z2dt * laplacu(2,:) * umask(2,:,1)316 IF( (nbondi == 1) .OR. (nbondi == 2)) spgu(nlci-2,:) = znugdt * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)317 IF( (nbondj == -1) .OR. (nbondj == 2) ) spgv(:,2) = znugdt * z2dt * laplacv(:,2) * vmask(:,2,1)318 IF( (nbondj == 1) .OR. (nbondj == 2)) spgv(:,nlcj-2) = znugdt * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)314 IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2 ,:) = znugdt * z2dt * laplacu(2 ,:) * umask(2 ,:,1) 315 IF( nbondi == 1 .OR. nbondi == 2 ) spgu(nlci-2,:) = znugdt * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1) 316 IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2 ) = znugdt * z2dt * laplacv(:,2 ) * vmask(: ,2,1) 317 IF( nbondj == 1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = znugdt * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1) 319 318 ENDIF 320 319 #endif … … 323 322 ! ( c a u t i o n : (ua,va) here are the after velocity not the 324 323 ! trend, the leap-frog time stepping will not 325 ! be done in dynnxt.F routine)324 ! be done in dynnxt.F90 routine) 326 325 DO jk = 1, jpkm1 327 326 DO jj = 2, jpjm1 … … 332 331 END DO 333 332 END DO 334 335 333 336 334 ! Sea surface elevation time stepping … … 358 356 ENDIF 359 357 360 ! ! print sum trends (used for debugging) 361 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg - ssh: ', mask1=tmask ) 362 358 ! write filtered free surface arrays in restart file 359 ! -------------------------------------------------- 360 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 361 362 ! print sum trends (used for debugging) 363 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg - ssh: ', mask1=tmask ) 364 ! 363 365 END SUBROUTINE dyn_spg_flt 366 367 368 SUBROUTINE flt_rst( kt, cdrw ) 369 !!--------------------------------------------------------------------- 370 !! *** ROUTINE ts_rst *** 371 !! 372 !! ** Purpose : Read or write filtered free surface arrays in restart file 373 !!---------------------------------------------------------------------- 374 INTEGER , INTENT(in) :: kt ! ocean time-step 375 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 376 !!---------------------------------------------------------------------- 377 378 IF( TRIM(cdrw) == 'READ' ) THEN 379 IF( iom_varid( numror, 'gcx' ) > 0 ) THEN 380 ! Caution : extra-hallow 381 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 382 CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) ) 383 CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) ) 384 CALL iom_get( numror, jpdom_local, 'sshb', sshb(:,:) ) 385 CALL iom_get( numror, jpdom_local, 'sshn', sshn(:,:) ) 386 IF( neuler == 0 ) THEN 387 sshb(:,:) = sshn(:,:) 388 gcxb(:,:) = gcx (:,:) 389 ENDIF 390 ELSE 391 gcx (:,:) = 0.e0 392 gcxb(:,:) = 0.e0 393 sshb(:,:) = 0.e0 394 sshn(:,:) = 0.e0 395 ENDIF 396 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 397 ! Caution : extra-hallow 398 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 399 CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) ) 400 CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) ) 401 CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:) ) 402 CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:) ) 403 ENDIF 404 ! 405 END SUBROUTINE flt_rst 364 406 365 407 #else
Note: See TracChangeset
for help on using the changeset viewer.
