Changeset 508 for trunk/NEMO/OPA_SRC/DYN
- Timestamp:
- 2006-10-03T17:58:55+02:00 (18 years ago)
- Location:
- trunk/NEMO/OPA_SRC/DYN
- Files:
-
- 4 edited
-
dynspg_flt.F90 (modified) (13 diffs)
-
dynspg_flt_jki.F90 (modified) (5 diffs)
-
dynspg_rl.F90 (modified) (15 diffs)
-
dynspg_ts.F90 (modified) (16 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 -
trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90
r503 r508 25 25 USE solfet ! FETI solver 26 26 USE solsor_e ! Successive Over-relaxation solver with MPP optimization 27 USE solver 27 28 USE obc_oce ! Lateral open boundary condition 28 29 USE obcdyn ! ocean open boundary condition (obc_dyn routines) … … 35 36 USE solmat ! matrix construction for elliptic solvers 36 37 USE agrif_opa_interp 38 USE restart ! only for lrst_oce 39 USE iom 37 40 38 41 IMPLICIT NONE … … 112 115 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~ (free surface constant volume, autotasking case)' 113 116 114 ! set to zero free surface specific arrays 115 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 116 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 117 ! set to zero free surface specific arrays 118 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 119 120 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 121 122 CALL solver_init( nit000 ) ! Elliptic solver initialisation 123 124 ! read filtered free surface arrays in restart file 125 CALL flt_rst( nit000, 'READ' ) ! read or initialize the following fields: 126 ! ! gcx, gcxb, sshb, sshn 117 127 ENDIF 118 128 … … 354 364 CALL lbc_lnk( sshn, 'T', 1. ) 355 365 366 ! write filtered free surface arrays in restart file 367 ! -------------------------------------------------- 368 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 369 356 370 IF(ln_ctl) THEN ! print sum trends (used for debugging) 357 371 CALL prt_ctl( tab3d_1=ua , clinfo1=' spg - Ua : ', mask1=umask, & … … 362 376 END SUBROUTINE dyn_spg_flt_jki 363 377 378 SUBROUTINE flt_rst( kt, cdrw ) 379 !!--------------------------------------------------------------------- 380 !! *** ROUTINE ts_rst *** 381 !! 382 !! ** Purpose : Read or write filtered free surface arrays in restart file 383 !!---------------------------------------------------------------------- 384 INTEGER , INTENT(in) :: kt ! ocean time-step 385 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 386 !!---------------------------------------------------------------------- 387 388 IF( TRIM(cdrw) == 'READ' ) THEN 389 IF( iom_varid( numror, 'gcx' ) > 0 ) THEN 390 ! Caution : extra-hallow 391 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 392 CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) ) 393 CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) ) 394 CALL iom_get( numror, jpdom_local, 'sshb', sshb(:,:) ) 395 CALL iom_get( numror, jpdom_local, 'sshn', sshn(:,:) ) 396 IF( neuler == 0 ) THEN 397 sshb(:,:) = sshn(:,:) 398 gcxb(:,:) = gcx (:,:) 399 ENDIF 400 ELSE 401 gcx (:,:) = 0.e0 402 gcxb(:,:) = 0.e0 403 sshb(:,:) = 0.e0 404 sshn(:,:) = 0.e0 405 ENDIF 406 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 407 ! Caution : extra-hallow 408 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 409 CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) ) 410 CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) ) 411 CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:) ) 412 CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:) ) 413 ENDIF 414 ! 415 END SUBROUTINE flt_rst 416 364 417 #else 365 418 !!---------------------------------------------------------------------- -
trunk/NEMO/OPA_SRC/DYN/dynspg_rl.F90
r474 r508 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 !! History : 7.0 ! 96-05 (G. Madec, M. Imbard, M. Guyon) rewitting in 1 7 !! routine, without pointers, and with the same matrix 8 !! for sor and pcg, mpp exchange, and symmetric conditions 9 !! " " ! 96-07 (A. Weaver) Euler forward step 10 !! " " ! 96-11 (A. Weaver) correction to preconditioning: 11 !! 8.0 ! 98-02 (M. Guyon) FETI method 12 !! " " ! 97-09 (J.-M. Molines) Open boundaries 13 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 14 !! ! 02-11 (C. Talandier, A-M Treguier) Open boundaries 15 !! 9.0 ! 04-08 (C. Talandier) New trends organization 16 !! " " ! 05-11 (V. Garnier) Surface pressure gradient organization 17 !! " " ! 06-07 (S. Masson) distributed restart using iom 18 !!--------------------------------------------------------------------- 6 19 #if defined key_dynspg_rl || defined key_esopa 7 20 !!---------------------------------------------------------------------- … … 10 23 !! dyn_spg_rl : update the momentum trend with the surface pressure 11 24 !! for the rigid-lid case. 12 !! ----------------------------------------------------------------------13 !! * Modules used25 !! rl_rst : read/write the rigid-lid restart fields in the ocean restart file 26 !!---------------------------------------------------------------------- 14 27 USE oce ! ocean dynamics and tracers 15 28 USE dom_oce ! ocean space and time domain … … 19 32 USE zdf_oce ! ocean vertical physics 20 33 USE sol_oce ! ocean elliptic solver 34 USE solver ! solver initialization 21 35 USE solpcg ! preconditionned conjugate gradient solver 22 36 USE solsor ! Successive Over-relaxation solver … … 28 42 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 43 USE in_out_manager ! I/O manager 44 USE iom 45 USE restart ! only for lrst_oce 30 46 31 47 IMPLICIT NONE 32 48 PRIVATE 33 49 34 !! * Accessibility35 50 PUBLIC dyn_spg_rl ! called by step.F90 36 51 … … 42 57 !! OPA 9.0 , LOCEAN-IPSL (2005) 43 58 !! $Header$ 44 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt59 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 45 60 !!---------------------------------------------------------------------- 46 61 … … 76 91 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 77 92 !! 78 !! References : 79 !! Madec et al. 1988, ocean modelling, issue 78, 1-6. 80 !! 81 !! History : 82 !! ! 96-05 (G. Madec, M. Imbard, M. Guyon) rewitting in 1 83 !! routine, without pointers, and with the same matrix 84 !! for sor and pcg, mpp exchange, and symmetric conditions 85 !! ! 96-07 (A. Weaver) Euler forward step 86 !! ! 96-11 (A. Weaver) correction to preconditioning: 87 !! ! 98-02 (M. Guyon) FETI method 88 !! ! 98-05 (G. Roullet) free surface 89 !! ! 97-09 (J.-M. Molines) Open boundaries 90 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 91 !! ! 02-11 (C. Talandier, A-M Treguier) Open boundaries 92 !! 9.0 ! 04-08 (C. Talandier) New trends organization 93 !! " ! 05-11 (V. Garnier) Surface pressure gradient organization 93 !! References : Madec et al. 1988, ocean modelling, issue 78, 1-6. 94 94 !!--------------------------------------------------------------------- 95 !! * Arguments96 95 INTEGER, INTENT( in ) :: kt ! ocean time-step index 97 INTEGER, INTENT( out ) :: kindic ! solver flag, take a negative value 98 ! ! when the solver doesnot converge 99 !! * Local declarations 96 INTEGER, INTENT( out ) :: kindic ! solver flag (<0 when the solver does not converge) 100 97 INTEGER :: ji, jj, jk ! dummy loop indices 101 98 REAL(wp) :: zbsfa, zgcx, z2dt … … 114 111 115 112 ! set to zero rigid-lid specific arrays 116 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 117 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 118 ENDIF 119 120 ! 0. Initializations: 121 ! ------------------- 122 # if defined key_obc 123 ! space index on boundary arrays 124 ib = 1 125 ibm = 2 126 ibm2 = 3 127 ! time index on boundary arrays 128 it = 1 129 itm = 2 130 itm2 = 3 131 # endif 132 113 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 114 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 115 116 CALL solver_init( nit000 ) ! Elliptic solver initialisation 117 118 ! read rigid lid arrays in restart file 119 CALL rl_rst( nit000, 'READ' ) ! read or initialize the following fields: 120 ! ! gcx, gcxb, bsfb, bsfn, bsfd 121 ENDIF 122 123 ! Vertically averaged momentum trend 124 ! ------------------------------------ 133 125 ! ! =============== 134 126 DO jj = 2, jpjm1 ! Vertical slab 135 127 ! ! =============== 136 137 ! 1. Vertically averaged momentum trend 138 ! ------------------------------------- 139 ! initialization to zero 140 spgu(:,jj) = 0. 128 129 spgu(:,jj) = 0. ! initialization to zero 141 130 spgv(:,jj) = 0. 142 ! vertical sum 143 DO jk = 1, jpk 131 DO jk = 1, jpk ! vertical sum 144 132 DO ji = 2, jpim1 145 133 spgu(ji,jj) = spgu(ji,jj) + ua(ji,jj,jk) * fse3u(ji,jj,jk) * umask(ji,jj,jk) … … 147 135 END DO 148 136 END DO 149 ! divide by the depth 150 spgu(:,jj) = spgu(:,jj) * hur(:,jj) 137 spgu(:,jj) = spgu(:,jj) * hur(:,jj) ! divide by the depth 151 138 spgv(:,jj) = spgv(:,jj) * hvr(:,jj) 152 139 … … 155 142 ! ! =============== 156 143 157 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,158 159 144 ! Boundary conditions on (spgu,spgv) 160 145 CALL lbc_lnk( spgu, 'U', -1. ) 161 146 CALL lbc_lnk( spgv, 'V', -1. ) 162 147 163 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 164 165 ! 2. Barotropic streamfunction trend (bsfd) 148 ! Barotropic streamfunction trend (bsfd) 166 149 ! ---------------------------------- 167 168 ! Curl of the vertically averaged velocity 169 DO jj = 2, jpjm1 150 DO jj = 2, jpjm1 ! Curl of the vertically averaged velocity 170 151 DO ji = 2, jpim1 171 152 gcb(ji,jj) = -gcdprc(ji,jj) & 172 173 153 & *( ( e2v(ji+1,jj )*spgv(ji+1,jj ) - e2v(ji,jj)*spgv(ji,jj) ) & 154 & -( e1u(ji ,jj+1)*spgu(ji ,jj+1) - e1u(ji,jj)*spgu(ji,jj) ) ) 174 155 END DO 175 156 END DO 176 157 177 158 # if defined key_obc 178 ! Open boundary contribution 179 DO jj = 2, jpjm1 159 DO jj = 2, jpjm1 ! Open boundary contribution 180 160 DO ji = 2, jpim1 181 161 gcb(ji,jj) = gcb(ji,jj) - gcdprc(ji,jj) * gcbob(ji,jj) … … 198 178 ! applied the lateral boundary conditions 199 179 IF( nsolv == 4) CALL lbc_lnk_e( gcb, c_solver_pt, 1. ) 200 201 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,202 180 203 181 ! Relative precision (computation on one processor) … … 234 212 ENDIF 235 213 236 237 214 ! bsf trend update 238 215 ! ---------------- 239 240 216 bsfd(1:nlci,1:nlcj) = gcx(1:nlci,1:nlcj) 241 242 217 243 218 ! update bsf trend with islands trend 244 219 ! ----------------------------------- 245 246 220 IF( lk_isl ) CALL isl_dyn_spg ! update bsfd 247 248 221 249 222 # if defined key_obc … … 337 310 ! 4. Barotropic stream function and array swap 338 311 ! -------------------------------------------- 339 340 312 ! Leap-frog time scheme, time filter and array swap 341 313 IF( neuler == 0 .AND. kt == nit000 ) THEN … … 362 334 363 335 # if defined key_obc 336 ib = 1 ! space index on boundary arrays 337 ibm = 2 338 ibm2 = 3 339 it = 1 ! time index on boundary arrays 340 itm = 2 341 itm2 = 3 342 364 343 ! Swap of boundary arrays 365 344 IF( lp_obc_east ) THEN … … 499 478 ENDIF 500 479 # endif 501 !502 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,503 480 504 481 ! add the surface pressure trend to the general trend 505 482 ! ----------------------------------------------------- 506 507 483 DO jj = 2, jpjm1 508 509 484 ! update the surface pressure gradient with the barotropic trend 510 485 DO ji = 2, jpim1 … … 519 494 END DO 520 495 END DO 521 522 END DO 496 END DO 497 498 ! write rigid lid arrays in restart file 499 ! -------------------------------------- 500 IF( lrst_oce ) CALL rl_rst( kt, 'WRITE' ) 523 501 524 502 END SUBROUTINE dyn_spg_rl 503 504 505 SUBROUTINE rl_rst( kt, cdrw ) 506 !!--------------------------------------------------------------------- 507 !! *** ROUTINE rl_rst *** 508 !! 509 !! ** Purpose : Read or write rigid-lid arrays in ocean restart file 510 !!---------------------------------------------------------------------- 511 INTEGER , INTENT(in) :: kt ! ocean time-step 512 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 513 !!---------------------------------------------------------------------- 514 ! 515 IF( TRIM(cdrw) == 'READ' ) THEN 516 IF( iom_varid( numror, 'gcx' ) > 0 ) THEN 517 ! Caution : extra-hallow 518 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 519 CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) ) 520 CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) ) 521 CALL iom_get( numror, jpdom_local, 'bsfb', bsfb(:,:) ) 522 CALL iom_get( numror, jpdom_local, 'bsfn', bsfn(:,:) ) 523 CALL iom_get( numror, jpdom_local, 'bsfd', bsfd(:,:) ) 524 IF( neuler == 0 ) THEN 525 gcxb(:,:) = gcx (:,:) 526 bsfb(:,:) = bsfn(:,:) 527 ENDIF 528 ELSE 529 gcx (:,:) = 0.e0 530 gcxb(:,:) = 0.e0 531 bsfb(:,:) = 0.e0 532 bsfn(:,:) = 0.e0 533 bsfd(:,:) = 0.e0 534 ENDIF 535 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 536 ! Caution : extra-hallow, gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 537 CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) ) 538 CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) ) 539 CALL iom_rstput( kt, nitrst, numrow, 'bsfb', bsfb(:,:) ) 540 CALL iom_rstput( kt, nitrst, numrow, 'bsfn', bsfn(:,:) ) 541 CALL iom_rstput( kt, nitrst, numrow, 'bsfd', bsfd(:,:) ) 542 ENDIF 543 ! 544 END SUBROUTINE rl_rst 525 545 526 546 #else -
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r455 r508 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 !! History : 9.0 ! 04-12 (L. Bessieres, G. Madec) Original code 7 !! " " ! 05-11 (V. Garnier, G. Madec) optimization 8 !! 9.0 ! 06-08 (S. Masson) distributed restart using iom 9 !!--------------------------------------------------------------------- 6 10 #if ( defined key_dynspg_ts && ! defined key_mpp_omp ) || defined key_esopa 7 11 !!---------------------------------------------------------------------- … … 9 13 !! NOT 'key_mpp_omp' k-j-i loop (vector opt.) 10 14 !!---------------------------------------------------------------------- 15 !!---------------------------------------------------------------------- 11 16 !! dyn_spg_ts : compute surface pressure gradient trend using a time- 12 17 !! splitting scheme and add to the general trend 18 !! ts_rst : read/write the time-splitting restart fields in the ocean restart file 13 19 !!---------------------------------------------------------------------- 14 20 !! * Modules used … … 27 33 USE dynspg_oce ! surface pressure gradient variables 28 34 USE in_out_manager ! I/O manager 35 USE iom 36 USE restart ! only for lrst_oce 29 37 30 38 IMPLICIT NONE 31 39 PRIVATE 32 40 33 !! * Accessibility34 41 PUBLIC dyn_spg_ts ! routine called by step.F90 42 43 REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne, & ! triad of coriolis parameter 44 & ftsw, ftse ! (only used with een vorticity scheme) 45 35 46 36 47 !! * Substitutions … … 74 85 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 75 86 !! 76 !! References : 77 !! Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 78 !! 79 !! History : 80 !! 9.0 ! 04-12 (L. Bessieres, G. Madec) Original code 81 !! ! 05-11 (V. Garnier, G. Madec) optimization 87 !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 82 88 !!--------------------------------------------------------------------- 83 !! * Arguments84 89 INTEGER, INTENT( in ) :: kt ! ocean time-step index 85 90 … … 97 102 zsshb_e, zub_e, zvb_e, & ! " " 98 103 zun_e, zvn_e ! " " 99 REAL(wp), DIMENSION(jpi,jpj),SAVE :: &100 ztnw, ztne, ztsw, ztse101 104 !!---------------------------------------------------------------------- 102 105 … … 109 112 110 113 IF( kt == nit000 ) THEN 111 114 ! 112 115 IF(lwp) WRITE(numout,*) 113 116 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend' … … 115 118 IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', FLOOR( 2*rdt/rdtbt ) 116 119 117 IF( .NOT. ln_rstart ) THEN 118 ! initialize barotropic specific arrays 119 sshb_b(:,:) = sshb(:,:) 120 sshn_b(:,:) = sshn(:,:) 121 un_b(:,:) = 0.e0 122 vn_b(:,:) = 0.e0 123 ! vertical sum 124 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 125 DO jk = 1, jpkm1 126 DO ji = 1, jpij 127 un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 128 vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 129 END DO 130 END DO 131 ELSE ! No vector opt. 132 DO jk = 1, jpkm1 133 un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) 134 vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 135 END DO 136 ENDIF 137 ENDIF 120 CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: 121 ! ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 122 138 123 ssha_e(:,:) = sshn(:,:) 139 124 ua_e(:,:) = un_b(:,:) … … 141 126 142 127 IF( ln_dynvor_een ) THEN 143 ztne(1,:) = 0.e0 ; ztnw(1,:) = 0.e0 ; ztse(1,:) = 0.e0 ; ztsw(1,:) = 0.e0128 ftne(1,:) = 0.e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0 144 129 DO jj = 2, jpj 145 130 DO ji = fs_2, jpi ! vector opt. 146 ztne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3.147 ztnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3.148 ztse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3.149 ztsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3.131 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3. 132 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3. 133 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3. 134 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3. 150 135 END DO 151 136 END DO 152 137 ENDIF 153 138 ! 154 139 ENDIF 155 140 156 141 ! Local constant initialization 157 142 ! -------------------------------- … … 216 201 END DO 217 202 END DO 218 203 ! 219 204 ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme 220 205 DO jj = 2, jpjm1 … … 228 213 END DO 229 214 END DO 230 215 ! 231 216 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 232 217 zfac25 = 0.25 … … 241 226 END DO 242 227 END DO 243 228 ! 244 229 ENDIF 245 230 … … 300 285 DO jit = 1, icycle ! sub-time-step loop ! 301 286 ! ! ==================== ! 302 303 287 z2dt_e = 2. * rdtbt 304 288 IF ( jit == 1 ) z2dt_e = rdtbt … … 360 344 END DO 361 345 END DO 362 346 ! 363 347 ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme 364 348 DO jj = 2, jpjm1 … … 379 363 END DO 380 364 END DO 381 365 ! 382 366 ELSEIF ( ln_dynvor_een ) THEN ! energy and enstrophy conserving scheme 383 367 zfac25 = 0.25 … … 397 381 END DO 398 382 END DO 383 ! 399 384 ENDIF 400 385 … … 504 489 END DO 505 490 506 IF(ln_ctl) THEN ! print sum trends (used for debugging) 507 CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask) 491 ! write filtered free surface arrays in restart file 492 ! -------------------------------------------------- 493 IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) 494 495 ! print sum trends (used for debugging) 496 IF( ln_ctl ) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask ) 497 ! 498 END SUBROUTINE dyn_spg_ts 499 500 501 SUBROUTINE ts_rst( kt, cdrw ) 502 !!--------------------------------------------------------------------- 503 !! *** ROUTINE ts_rst *** 504 !! 505 !! ** Purpose : Read or write time-splitting arrays in restart file 506 !!---------------------------------------------------------------------- 507 INTEGER , INTENT(in) :: kt ! ocean time-step 508 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 509 ! 510 INTEGER :: ji, jk ! dummy loop indices 511 !!---------------------------------------------------------------------- 512 ! 513 IF( TRIM(cdrw) == 'READ' ) THEN 514 IF( iom_varid( numror, 'sshn' ) > 0 ) THEN 515 CALL iom_get( numror, jpdom_local, 'sshb' , sshb(:,:) ) 516 CALL iom_get( numror, jpdom_local, 'sshn' , sshn(:,:) ) 517 IF( neuler == 0 ) sshb(:,:) = sshn(:,:) 518 ELSE 519 sshb(:,:) = 0.e0 520 sshn(:,:) = 0.e0 521 ENDIF 522 IF( iom_varid( numror, 'sshn_b' ) > 0 ) THEN 523 CALL iom_get( numror, jpdom_local, 'sshb_b', sshb_b(:,:) ) ! free surface issued 524 CALL iom_get( numror, jpdom_local, 'sshn_b', sshn_b(:,:) ) ! from time-splitting loop 525 CALL iom_get( numror, jpdom_local, 'un_b' , un_b (:,:) ) ! horizontal transports issued 526 CALL iom_get( numror, jpdom_local, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 527 IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:) 528 ELSE 529 sshb_b(:,:) = sshb(:,:) 530 sshn_b(:,:) = sshn(:,:) 531 un_b (:,:) = 0.e0 532 vn_b (:,:) = 0.e0 533 ! vertical sum 534 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 535 DO jk = 1, jpkm1 536 DO ji = 1, jpij 537 un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 538 vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 539 END DO 540 END DO 541 ELSE ! No vector opt. 542 DO jk = 1, jpkm1 543 un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) 544 vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 545 END DO 546 ENDIF 547 ENDIF 548 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 549 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb (:,:) ) 550 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn (:,:) ) 551 CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) ) ! free surface issued 552 CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) ) ! from barotropic loop 553 CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! horizontal transports issued 554 CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 508 555 ENDIF 509 510 END SUBROUTINE dyn_spg_ts 556 ! 557 END SUBROUTINE ts_rst 558 511 559 #else 512 560 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
