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 branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 4328

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

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

Download in other formats: