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/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago
#1260 : see the associated wiki page for explanation
Property svn:keywords set to `Id`
File size: 18.0 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.
[4616]	168	spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) * r1_e1u(ji,jj)
	169	spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) * r1_e2v(ji,jj)
[592]	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	! ---------------------------------------------
[4616]	194	! vertical sum
[358]	195	DO jj = 2, jpjm1
	196	DO ji = fs_2, fs_jpim1 ! vector opt.
[4616]	197	spgu(ji,jj) = fse3u_a(ji,jj,1) * ua(ji,jj,1)
	198	spgv(ji,jj) = fse3v_a(ji,jj,1) * va(ji,jj,1)
[358]	199	END DO
	200	END DO
[4616]	201	DO jk = 2, jpkm1
	202	DO jj = 2, jpjm1
	203	DO ji = 2, jpim1
	204	spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk)
	205	spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk)
[358]	206	END DO
	207	END DO
[4616]	208	END DO
[358]	209
	210	! transport: multiplied by the horizontal scale factor
	211	DO jj = 2, jpjm1
	212	DO ji = fs_2, fs_jpim1 ! vector opt.
	213	spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
	214	spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
	215	END DO
	216	END DO
[1601]	217	CALL lbc_lnk( spgu, 'U', -1. ) ! lateral boundary conditions
[358]	218	CALL lbc_lnk( spgv, 'V', -1. )
	219
[1438]	220	IF( lk_vvl ) CALL sol_mat( kt ) ! build the matrix at kt (vvl case only)
[592]	221
[358]	222	! Right hand side of the elliptic equation and first guess
[1601]	223	! --------------------------------------------------------
[358]	224	DO jj = 2, jpjm1
	225	DO ji = fs_2, fs_jpim1 ! vector opt.
	226	! Divergence of the after vertically averaged velocity
	227	zgcb = spgu(ji,jj) - spgu(ji-1,jj) &
	228	+ spgv(ji,jj) - spgv(ji,jj-1)
	229	gcb(ji,jj) = gcdprc(ji,jj) * zgcb
	230	! First guess of the after barotropic transport divergence
	231	zbtd = gcx(ji,jj)
	232	gcx (ji,jj) = 2. * zbtd - gcxb(ji,jj)
	233	gcxb(ji,jj) = zbtd
	234	END DO
	235	END DO
	236	! applied the lateral boundary conditions
[3609]	237	IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 ) CALL lbc_lnk_e( gcb, c_solver_pt, 1., jpr2di, jpr2dj )
[358]	238
[392]	239	#if defined key_agrif
[413]	240	IF( .NOT. AGRIF_ROOT() ) THEN
[389]	241	! add contribution of gradient of after barotropic transport divergence
[508]	242	IF( nbondi == -1 .OR. nbondi == 2 ) gcb(3 ,:) = &
[1601]	243	& gcb(3 ,:) - z2dtg * z2dt * laplacu(2 ,:) * gcdprc(3 ,:) * hu(2 ,:) * e2u(2 ,:)
[508]	244	IF( nbondi == 1 .OR. nbondi == 2 ) gcb(nlci-2,:) = &
[1601]	245	& gcb(nlci-2,:) + z2dtg * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)
[508]	246	IF( nbondj == -1 .OR. nbondj == 2 ) gcb(: ,3) = &
[1601]	247	& gcb(:,3 ) - z2dtg * z2dt * laplacv(:,2 ) * gcdprc(:,3 ) * hv(:,2 ) * e1v(:,2 )
[508]	248	IF( nbondj == 1 .OR. nbondj == 2 ) gcb(:,nlcj-2) = &
[1601]	249	& gcb(:,nlcj-2) + z2dtg * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)
[413]	250	ENDIF
[389]	251	#endif
	252
	253
[358]	254	! Relative precision (computation on one processor)
	255	! ------------------
[1438]	256	rnorme =0.e0
[2528]	257	rnorme = GLOB_SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) )
[358]	258
	259	epsr = eps * eps * rnorme
	260	ncut = 0
[508]	261	! if rnorme is 0, the solution is 0, the solver is not called
[2715]	262	IF( rnorme == 0._wp ) THEN
	263	gcx(:,:) = 0._wp
	264	res = 0._wp
[358]	265	niter = 0
	266	ncut = 999
	267	ENDIF
	268
	269	! Evaluate the next transport divergence
	270	! --------------------------------------
	271	! Iterarive solver for the elliptic equation (except IF sol.=0)
	272	! (output in gcx with boundary conditions applied)
	273	kindic = 0
	274	IF( ncut == 0 ) THEN
[1601]	275	IF ( nn_solv == 1 ) THEN ; CALL sol_pcg( kindic ) ! diagonal preconditioned conjuguate gradient
	276	ELSEIF( nn_solv == 2 ) THEN ; CALL sol_sor( kindic ) ! successive-over-relaxation
[358]	277	ENDIF
	278	ENDIF
	279
	280	! Transport divergence gradient multiplied by z2dt
	281	! --------------------------------------------====
	282	DO jj = 2, jpjm1
	283	DO ji = fs_2, fs_jpim1 ! vector opt.
	284	! trend of Transport divergence gradient
[4616]	285	ztdgu = z2dtg * (gcx(ji+1,jj ) - gcx(ji,jj) ) * r1_e1u(ji,jj)
	286	ztdgv = z2dtg * (gcx(ji ,jj+1) - gcx(ji,jj) ) * r1_e2v(ji,jj)
[358]	287	! multiplied by z2dt
[4328]	288	#if defined key_bdy
[3765]	289	IF(lk_bdy) THEN
[911]	290	! caution : grad D = 0 along open boundaries
[3765]	291	spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
	292	spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)
	293	ELSE
	294	spgu(ji,jj) = z2dt * ztdgu
	295	spgv(ji,jj) = z2dt * ztdgv
	296	ENDIF
[358]	297	#else
	298	spgu(ji,jj) = z2dt * ztdgu
	299	spgv(ji,jj) = z2dt * ztdgv
	300	#endif
	301	END DO
	302	END DO
	303
[1876]	304	#if defined key_agrif
[413]	305	IF( .NOT. Agrif_Root() ) THEN
	306	! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
[1601]	307	IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2 ,:) = z2dtg * z2dt * laplacu(2 ,:) * umask(2 ,:,1)
	308	IF( nbondi == 1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z2dtg * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
	309	IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2 ) = z2dtg * z2dt * laplacv(:,2 ) * vmask(: ,2,1)
	310	IF( nbondj == 1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z2dtg * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
[389]	311	ENDIF
[1876]	312	#endif
[1438]	313	! Add the trends multiplied by z2dt to the after velocity
	314	! -------------------------------------------------------
[358]	315	! ( c a u t i o n : (ua,va) here are the after velocity not the
	316	! trend, the leap-frog time stepping will not
[508]	317	! be done in dynnxt.F90 routine)
[358]	318	DO jk = 1, jpkm1
	319	DO jj = 2, jpjm1
	320	DO ji = fs_2, fs_jpim1 ! vector opt.
[1438]	321	ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk)
	322	va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk)
[358]	323	END DO
	324	END DO
	325	END DO
	326
[508]	327	! write filtered free surface arrays in restart file
	328	! --------------------------------------------------
	329	IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
	330	!
[3294]	331	!
	332	IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_flt')
	333	!
[358]	334	END SUBROUTINE dyn_spg_flt
	335
[508]	336
	337	SUBROUTINE flt_rst( kt, cdrw )
[2715]	338	!!---------------------------------------------------------------------
	339	!! * ROUTINE ts_rst *
	340	!!
	341	!! ** Purpose : Read or write filtered free surface arrays in restart file
	342	!!----------------------------------------------------------------------
	343	INTEGER , INTENT(in) :: kt ! ocean time-step
	344	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	345	!!----------------------------------------------------------------------
	346	!
	347	IF( TRIM(cdrw) == 'READ' ) THEN
	348	IF( iom_varid( numror, 'gcx', ldstop = .FALSE. ) > 0 ) THEN
[508]	349	! Caution : extra-hallow
	350	! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
[2715]	351	CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
	352	CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
	353	IF( neuler == 0 ) gcxb(:,:) = gcx (:,:)
	354	ELSE
	355	gcx (:,:) = 0.e0
	356	gcxb(:,:) = 0.e0
	357	ENDIF
	358	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
[508]	359	! Caution : extra-hallow
	360	! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
[2715]	361	CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
	362	CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
	363	ENDIF
	364	!
[508]	365	END SUBROUTINE flt_rst
	366
[358]	367	#else
	368	!!----------------------------------------------------------------------
	369	!! Default case : Empty module No standart free surface cst volume
	370	!!----------------------------------------------------------------------
	371	CONTAINS
	372	SUBROUTINE dyn_spg_flt( kt, kindic ) ! Empty routine
	373	WRITE(,) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
	374	END SUBROUTINE dyn_spg_flt
[657]	375	SUBROUTINE flt_rst ( kt, cdrw ) ! Empty routine
	376	INTEGER , INTENT(in) :: kt ! ocean time-step
	377	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	378	WRITE(,) 'flt_rst: You should not have seen this print! error?', kt, cdrw
	379	END SUBROUTINE flt_rst
[358]	380	#endif
	381
	382	!!======================================================================
	383	END MODULE dynspg_flt

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

Download in other formats: