New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_flt.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 3260

Last change on this file since 3260 was 2715, checked in by rblod, 13 years ago
First attempt to put dynamic allocation on the trunk
Property svn:keywords set to `Id`
File size: 18.9 KB

Rev	Line
[358]	1	MODULE dynspg_flt
	2	!!======================================================================
	3	!! * MODULE dynspg_flt *
	4	!! Ocean dynamics: surface pressure gradient trend
	5	!!======================================================================
[1438]	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
[358]	15	!!----------------------------------------------------------------------
[575]	16	#if defined key_dynspg_flt \|\| defined key_esopa
[508]	17	!!----------------------------------------------------------------------
[358]	18	!! 'key_dynspg_flt' filtered free surface
	19	!!----------------------------------------------------------------------
[1601]	20	!! dyn_spg_flt : update the momentum trend with the surface pressure gradient in the filtered free surface case
[508]	21	!! flt_rst : read/write the time-splitting restart fields in the ocean restart file
[358]	22	!!----------------------------------------------------------------------
	23	USE oce ! ocean dynamics and tracers
	24	USE dom_oce ! ocean space and time domain
	25	USE zdf_oce ! ocean vertical physics
[888]	26	USE sbc_oce ! surface boundary condition: ocean
	27	USE obc_oce ! Lateral open boundary condition
	28	USE sol_oce ! ocean elliptic solver
[719]	29	USE phycst ! physical constants
[888]	30	USE domvvl ! variable volume
[1683]	31	USE dynadv ! advection
[1601]	32	USE solmat ! matrix construction for elliptic solvers
[358]	33	USE solpcg ! preconditionned conjugate gradient solver
	34	USE solsor ! Successive Over-relaxation solver
	35	USE obcdyn ! ocean open boundary condition (obc_dyn routines)
	36	USE obcvol ! ocean open boundary condition (obc_vol routines)
[911]	37	USE bdy_oce ! Unstructured open boundaries condition
	38	USE bdydyn ! Unstructured open boundaries condition (bdy_dyn routine)
	39	USE bdyvol ! Unstructured open boundaries condition (bdy_vol routine)
[2528]	40	USE cla ! cross land advection
[888]	41	USE in_out_manager ! I/O manager
[358]	42	USE lib_mpp ! distributed memory computing library
	43	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	44	USE prtctl ! Print control
[508]	45	USE iom
	46	USE restart ! only for lrst_oce
[2528]	47	USE lib_fortran
	48	#if defined key_agrif
	49	USE agrif_opa_interp
	50	#endif
[358]	51
	52	IMPLICIT NONE
	53	PRIVATE
	54
[1438]	55	PUBLIC dyn_spg_flt ! routine called by step.F90
	56	PUBLIC flt_rst ! routine called by istate.F90
[358]	57
	58	!! * Substitutions
	59	# include "domzgr_substitute.h90"
	60	# include "vectopt_loop_substitute.h90"
	61	!!----------------------------------------------------------------------
[2528]	62	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	63	!! $Id$
[2715]	64	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[358]	65	!!----------------------------------------------------------------------
	66	CONTAINS
	67
	68	SUBROUTINE dyn_spg_flt( kt, kindic )
	69	!!----------------------------------------------------------------------
	70	!! * routine dyn_spg_flt *
	71	!!
	72	!! ** Purpose : Compute the now trend due to the surface pressure
	73	!! gradient in case of filtered free surface formulation and add
	74	!! it to the general trend of momentum equation.
	75	!!
	76	!! ** Method : Filtered free surface formulation. The surface
	77	!! pressure gradient is given by:
[1601]	78	!! spgu = 1/rau0 d/dx(ps) = 1/e1u di( sshn + btda )
	79	!! spgv = 1/rau0 d/dy(ps) = 1/e2v dj( sshn + btda )
[358]	80	!! where sshn is the free surface elevation and btda is the after
[1438]	81	!! time derivative of the free surface elevation
	82	!! -1- evaluate the surface presure trend (including the addi-
[358]	83	!! tional force) in three steps:
	84	!! a- compute the right hand side of the elliptic equation:
	85	!! gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
	86	!! where (spgu,spgv) are given by:
	87	!! spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
[2528]	88	!! - grav 2 rdt hu /e1u di[sshn + (emp-rnf)]
[358]	89	!! spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
[2528]	90	!! - grav 2 rdt hv /e2v dj[sshn + (emp-rnf)]
[358]	91	!! and define the first guess from previous computation :
	92	!! zbtd = btda
	93	!! btda = 2 zbtd - btdb
	94	!! btdb = zbtd
	95	!! b- compute the relative accuracy to be reached by the
	96	!! iterative solver
	97	!! c- apply the solver by a call to sol... routine
[1438]	98	!! -2- compute and add the free surface pressure gradient inclu-
[358]	99	!! ding the additional force used to stabilize the equation.
	100	!!
	101	!! ** Action : - Update (ua,va) with the surf. pressure gradient trend
	102	!!
[508]	103	!! References : Roullet and Madec 1999, JGR.
[358]	104	!!---------------------------------------------------------------------
[2715]	105	USE oce, ONLY: zub => ta , zvb => sa ! (ta,sa) used as workspace
[1438]	106	!!
[1601]	107	INTEGER, INTENT(in ) :: kt ! ocean time-step index
	108	INTEGER, INTENT( out) :: kindic ! solver convergence flag (<0 if not converge)
[508]	109	!!
[2715]	110	INTEGER :: ji, jj, jk ! dummy loop indices
	111	REAL(wp) :: z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv ! local scalars
[358]	112	!!----------------------------------------------------------------------
[508]	113	!
[358]	114	IF( kt == nit000 ) THEN
	115	IF(lwp) WRITE(numout,*)
	116	IF(lwp) WRITE(numout,*) 'dyn_spg_flt : surface pressure gradient trend'
	117	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ (free surface constant volume case)'
	118
	119	! set to zero free surface specific arrays
[2715]	120	spgu(:,:) = 0._wp ! surface pressure gradient (i-direction)
	121	spgv(:,:) = 0._wp ! surface pressure gradient (j-direction)
[508]	122
	123	! read filtered free surface arrays in restart file
[1200]	124	! when using agrif, sshn, gcx have to be read in istate
[1601]	125	IF(.NOT. lk_agrif) CALL flt_rst( nit000, 'READ' ) ! read or initialize the following fields:
[1438]	126	! ! gcx, gcxb
[358]	127	ENDIF
	128
	129	! Local constant initialization
[1601]	130	z2dt = 2. * rdt ! time step: leap-frog
	131	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest
	132	IF( neuler == 0 .AND. kt == nit000+1 ) CALL sol_mat( kt )
[358]	133	z2dtg = grav * z2dt
	134
[1683]	135	! Evaluate the masked next velocity (effect of the additional force not included)
	136	! ---------------------------------
	137	IF( lk_vvl ) THEN ! variable volume (surface pressure gradient already included in dyn_hpg)
	138	!
	139	IF( ln_dynadv_vec ) THEN ! vector form : applied on velocity
	140	DO jk = 1, jpkm1
	141	DO jj = 2, jpjm1
	142	DO ji = fs_2, fs_jpim1 ! vector opt.
	143	ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
	144	va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
	145	END DO
[592]	146	END DO
	147	END DO
[1683]	148	!
	149	ELSE ! flux form : applied on thickness weighted velocity
	150	DO jk = 1, jpkm1
	151	DO jj = 2, jpjm1
	152	DO ji = fs_2, fs_jpim1 ! vector opt.
	153	ua(ji,jj,jk) = ( ub(ji,jj,jk) * fse3u_b(ji,jj,jk) &
	154	& + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk) ) &
	155	& / fse3u_a(ji,jj,jk) * umask(ji,jj,jk)
	156	va(ji,jj,jk) = ( vb(ji,jj,jk) * fse3v_b(ji,jj,jk) &
	157	& + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk) ) &
	158	& / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk)
	159	END DO
	160	END DO
	161	END DO
	162	!
	163	ENDIF
	164	!
	165	ELSE ! fixed volume (add the surface pressure gradient + unweighted time stepping)
	166	!
	167	DO jj = 2, jpjm1 ! Surface pressure gradient (now)
[358]	168	DO ji = fs_2, fs_jpim1 ! vector opt.
[592]	169	spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
	170	spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
	171	END DO
	172	END DO
[1683]	173	DO jk = 1, jpkm1 ! unweighted time stepping
[592]	174	DO jj = 2, jpjm1
	175	DO ji = fs_2, fs_jpim1 ! vector opt.
[1438]	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)
[592]	178	END DO
[358]	179	END DO
	180	END DO
[1438]	181	!
[592]	182	ENDIF
	183
[358]	184	#if defined key_obc
[2528]	185	CALL obc_dyn( kt ) ! Update velocities on each open boundary with the radiation algorithm
	186	CALL obc_vol( kt ) ! Correction of the barotropic componant velocity to control the volume of the system
[358]	187	#endif
[911]	188	#if defined key_bdy
[2528]	189	CALL bdy_dyn_frs( kt ) ! Update velocities on unstructured boundary using the Flow Relaxation Scheme
	190	CALL bdy_vol( kt ) ! Correction of the barotropic component velocity to control the volume of the system
[911]	191	#endif
[392]	192	#if defined key_agrif
[508]	193	CALL Agrif_dyn( kt ) ! Update velocities on each coarse/fine interfaces
[389]	194	#endif
[2528]	195	IF( nn_cla == 1 ) CALL cla_dynspg( kt ) ! Cross Land Advection (update (ua,va))
[358]	196
	197	! compute the next vertically averaged velocity (effect of the additional force not included)
	198	! ---------------------------------------------
	199	DO jj = 2, jpjm1
	200	DO ji = fs_2, fs_jpim1 ! vector opt.
[2715]	201	spgu(ji,jj) = 0._wp
	202	spgv(ji,jj) = 0._wp
[358]	203	END DO
	204	END DO
	205
	206	! vertical sum
	207	!CDIR NOLOOPCHG
	208	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
	209	DO jk = 1, jpkm1
	210	DO ji = 1, jpij
	211	spgu(ji,1) = spgu(ji,1) + fse3u(ji,1,jk) * ua(ji,1,jk)
	212	spgv(ji,1) = spgv(ji,1) + fse3v(ji,1,jk) * va(ji,1,jk)
	213	END DO
	214	END DO
	215	ELSE ! No vector opt.
	216	DO jk = 1, jpkm1
	217	DO jj = 2, jpjm1
	218	DO ji = 2, jpim1
	219	spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)
	220	spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)
	221	END DO
	222	END DO
	223	END DO
	224	ENDIF
	225
	226	! transport: multiplied by the horizontal scale factor
	227	DO jj = 2, jpjm1
	228	DO ji = fs_2, fs_jpim1 ! vector opt.
	229	spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
	230	spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
	231	END DO
	232	END DO
[1601]	233	CALL lbc_lnk( spgu, 'U', -1. ) ! lateral boundary conditions
[358]	234	CALL lbc_lnk( spgv, 'V', -1. )
	235
[1438]	236	IF( lk_vvl ) CALL sol_mat( kt ) ! build the matrix at kt (vvl case only)
[592]	237
[358]	238	! Right hand side of the elliptic equation and first guess
[1601]	239	! --------------------------------------------------------
[358]	240	DO jj = 2, jpjm1
	241	DO ji = fs_2, fs_jpim1 ! vector opt.
	242	! Divergence of the after vertically averaged velocity
	243	zgcb = spgu(ji,jj) - spgu(ji-1,jj) &
	244	+ spgv(ji,jj) - spgv(ji,jj-1)
	245	gcb(ji,jj) = gcdprc(ji,jj) * zgcb
	246	! First guess of the after barotropic transport divergence
	247	zbtd = gcx(ji,jj)
	248	gcx (ji,jj) = 2. * zbtd - gcxb(ji,jj)
	249	gcxb(ji,jj) = zbtd
	250	END DO
	251	END DO
	252	! applied the lateral boundary conditions
[1601]	253	IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 ) CALL lbc_lnk_e( gcb, c_solver_pt, 1. )
[358]	254
[392]	255	#if defined key_agrif
[413]	256	IF( .NOT. AGRIF_ROOT() ) THEN
[389]	257	! add contribution of gradient of after barotropic transport divergence
[508]	258	IF( nbondi == -1 .OR. nbondi == 2 ) gcb(3 ,:) = &
[1601]	259	& gcb(3 ,:) - z2dtg * z2dt * laplacu(2 ,:) * gcdprc(3 ,:) * hu(2 ,:) * e2u(2 ,:)
[508]	260	IF( nbondi == 1 .OR. nbondi == 2 ) gcb(nlci-2,:) = &
[1601]	261	& gcb(nlci-2,:) + z2dtg * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)
[508]	262	IF( nbondj == -1 .OR. nbondj == 2 ) gcb(: ,3) = &
[1601]	263	& gcb(:,3 ) - z2dtg * z2dt * laplacv(:,2 ) * gcdprc(:,3 ) * hv(:,2 ) * e1v(:,2 )
[508]	264	IF( nbondj == 1 .OR. nbondj == 2 ) gcb(:,nlcj-2) = &
[1601]	265	& gcb(:,nlcj-2) + z2dtg * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)
[413]	266	ENDIF
[389]	267	#endif
	268
	269
[358]	270	! Relative precision (computation on one processor)
	271	! ------------------
[1438]	272	rnorme =0.e0
[2528]	273	rnorme = GLOB_SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) )
[358]	274
	275	epsr = eps * eps * rnorme
	276	ncut = 0
[508]	277	! if rnorme is 0, the solution is 0, the solver is not called
[2715]	278	IF( rnorme == 0._wp ) THEN
	279	gcx(:,:) = 0._wp
	280	res = 0._wp
[358]	281	niter = 0
	282	ncut = 999
	283	ENDIF
	284
	285	! Evaluate the next transport divergence
	286	! --------------------------------------
	287	! Iterarive solver for the elliptic equation (except IF sol.=0)
	288	! (output in gcx with boundary conditions applied)
	289	kindic = 0
	290	IF( ncut == 0 ) THEN
[1601]	291	IF ( nn_solv == 1 ) THEN ; CALL sol_pcg( kindic ) ! diagonal preconditioned conjuguate gradient
	292	ELSEIF( nn_solv == 2 ) THEN ; CALL sol_sor( kindic ) ! successive-over-relaxation
[358]	293	ENDIF
	294	ENDIF
	295
	296	! Transport divergence gradient multiplied by z2dt
	297	! --------------------------------------------====
	298	DO jj = 2, jpjm1
	299	DO ji = fs_2, fs_jpim1 ! vector opt.
	300	! trend of Transport divergence gradient
[1601]	301	ztdgu = z2dtg * (gcx(ji+1,jj ) - gcx(ji,jj) ) / e1u(ji,jj)
	302	ztdgv = z2dtg * (gcx(ji ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
[358]	303	! multiplied by z2dt
	304	#if defined key_obc
	305	! caution : grad D = 0 along open boundaries
[2528]	306	spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
	307	spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
[911]	308	#elif defined key_bdy
	309	! caution : grad D = 0 along open boundaries
	310	! Remark: The filtering force could be reduced here in the FRS zone
	311	! by multiplying spgu/spgv by (1-alpha) ??
	312	spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
	313	spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)
[358]	314	#else
	315	spgu(ji,jj) = z2dt * ztdgu
	316	spgv(ji,jj) = z2dt * ztdgv
	317	#endif
	318	END DO
	319	END DO
	320
[1876]	321	#if defined key_agrif
[413]	322	IF( .NOT. Agrif_Root() ) THEN
	323	! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
[1601]	324	IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2 ,:) = z2dtg * z2dt * laplacu(2 ,:) * umask(2 ,:,1)
	325	IF( nbondi == 1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z2dtg * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
	326	IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2 ) = z2dtg * z2dt * laplacv(:,2 ) * vmask(: ,2,1)
	327	IF( nbondj == 1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z2dtg * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
[389]	328	ENDIF
[1876]	329	#endif
[1438]	330	! Add the trends multiplied by z2dt to the after velocity
	331	! -------------------------------------------------------
[358]	332	! ( c a u t i o n : (ua,va) here are the after velocity not the
	333	! trend, the leap-frog time stepping will not
[508]	334	! be done in dynnxt.F90 routine)
[358]	335	DO jk = 1, jpkm1
	336	DO jj = 2, jpjm1
	337	DO ji = fs_2, fs_jpim1 ! vector opt.
[1438]	338	ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk)
	339	va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk)
[358]	340	END DO
	341	END DO
	342	END DO
	343
[508]	344	! write filtered free surface arrays in restart file
	345	! --------------------------------------------------
	346	IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
	347	!
[358]	348	END SUBROUTINE dyn_spg_flt
	349
[508]	350
	351	SUBROUTINE flt_rst( kt, cdrw )
[2715]	352	!!---------------------------------------------------------------------
	353	!! * ROUTINE ts_rst *
	354	!!
	355	!! ** Purpose : Read or write filtered free surface arrays in restart file
	356	!!----------------------------------------------------------------------
	357	INTEGER , INTENT(in) :: kt ! ocean time-step
	358	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	359	!!----------------------------------------------------------------------
	360	!
	361	IF( TRIM(cdrw) == 'READ' ) THEN
	362	IF( iom_varid( numror, 'gcx', ldstop = .FALSE. ) > 0 ) THEN
[508]	363	! Caution : extra-hallow
	364	! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
[2715]	365	CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
	366	CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
	367	IF( neuler == 0 ) gcxb(:,:) = gcx (:,:)
	368	ELSE
	369	gcx (:,:) = 0.e0
	370	gcxb(:,:) = 0.e0
	371	ENDIF
	372	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
[508]	373	! Caution : extra-hallow
	374	! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
[2715]	375	CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
	376	CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
	377	ENDIF
	378	!
[508]	379	END SUBROUTINE flt_rst
	380
[358]	381	#else
	382	!!----------------------------------------------------------------------
	383	!! Default case : Empty module No standart free surface cst volume
	384	!!----------------------------------------------------------------------
	385	CONTAINS
	386	SUBROUTINE dyn_spg_flt( kt, kindic ) ! Empty routine
	387	WRITE(,) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
	388	END SUBROUTINE dyn_spg_flt
[657]	389	SUBROUTINE flt_rst ( kt, cdrw ) ! Empty routine
	390	INTEGER , INTENT(in) :: kt ! ocean time-step
	391	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	392	WRITE(,) 'flt_rst: You should not have seen this print! error?', kt, cdrw
	393	END SUBROUTINE flt_rst
[358]	394	#endif
	395
	396	!!======================================================================
	397	END MODULE dynspg_flt

Note: See TracBrowser for help on using the repository browser.

Download in other formats: