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_rl.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/DYN/dynspg_rl.F90 @ 1191

Last change on this file since 1191 was 1152, checked in by rblod, 16 years ago
Convert cvs header to svn Id, step II
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 22.4 KB

Rev	Line
[3]	1	MODULE dynspg_rl
	2	!!======================================================================
	3	!! * MODULE dynspg_rl *
	4	!! Ocean dynamics: surface pressure gradient trend
	5	!!======================================================================
[508]	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	!!---------------------------------------------------------------------
[3]	19	#if defined key_dynspg_rl \|\| defined key_esopa
	20	!!----------------------------------------------------------------------
	21	!! 'key_dynspg_rl' rigid lid
	22	!!----------------------------------------------------------------------
	23	!! dyn_spg_rl : update the momentum trend with the surface pressure
	24	!! for the rigid-lid case.
[508]	25	!! rl_rst : read/write the rigid-lid restart fields in the ocean restart file
[3]	26	!!----------------------------------------------------------------------
	27	USE oce ! ocean dynamics and tracers
	28	USE dom_oce ! ocean space and time domain
	29	USE phycst ! physical constants
	30	USE ldftra_oce ! ocean active tracers: lateral physics
	31	USE ldfdyn_oce ! ocean dynamics: lateral physics
	32	USE zdf_oce ! ocean vertical physics
	33	USE sol_oce ! ocean elliptic solver
[508]	34	USE solver ! solver initialization
[3]	35	USE solpcg ! preconditionned conjugate gradient solver
	36	USE solsor ! Successive Over-relaxation solver
	37	USE solfet ! FETI solver
	38	USE solisl ! ???
	39	USE obc_oce ! Lateral open boundary condition
[216]	40	USE lib_mpp ! distributed memory computing library
	41	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[359]	42	USE in_out_manager ! I/O manager
[508]	43	USE iom
	44	USE restart ! only for lrst_oce
[3]	45
	46	IMPLICIT NONE
	47	PRIVATE
	48
	49	PUBLIC dyn_spg_rl ! called by step.F90
	50
	51	!! * Substitutions
	52	# include "domzgr_substitute.h90"
	53	# include "vectopt_loop_substitute.h90"
[31]	54	# include "obc_vectopt_loop_substitute.h90"
[3]	55	!!----------------------------------------------------------------------
[247]	56	!! OPA 9.0 , LOCEAN-IPSL (2005)
[1152]	57	!! $Id$
[508]	58	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[3]	59	!!----------------------------------------------------------------------
	60
	61	CONTAINS
	62
	63	SUBROUTINE dyn_spg_rl( kt, kindic )
	64	!!----------------------------------------------------------------------
	65	!! * routine dyn_spg_rl *
	66	!!
	67	!! ** Purpose : Compute the now trend due to the surface pressure
	68	!! gradient for the rigid-lid case, add it to the general trend of
	69	!! momentum equation.
	70	!!
	71	!! ** Method : Rigid-lid appromimation: the surface pressure gradient
	72	!! is given by:
	73	!! spgu = 1/rau0 d/dx(ps) = Mu + 1/(hu e2u) dj-1(bsfd)
	74	!! spgv = 1/rau0 d/dy(ps) = Mv - 1/(hv e1v) di-1(bsfd)
	75	!! where (Mu,Mv) is the vertically averaged momentum trend (i.e.
	76	!! the vertical ponderated sum of the general momentum trend),
	77	!! and bsfd is the barotropic streamfunction trend.
	78	!! The trend is computed as follows:
	79	!! -1- compute the vertically averaged momentum trend (Mu,Mv)
	80	!! -2- compute the barotropic streamfunction trend by solving an
	81	!! ellipic equation using a diagonal preconditioned conjugate
	82	!! gradient or a successive-over-relaxation method (depending
	83	!! on nsolv, a namelist parameter).
[78]	84	!! -3- add to bsfd the island trends if lk_isl=T
[3]	85	!! -4- compute the after streamfunction is for further diagnos-
	86	!! tics using a leap-frog scheme.
	87	!! -5- add the momentum trend associated with the surface pres-
	88	!! sure gradient to the general trend.
	89	!!
	90	!! ** Action : - Update (ua,va) with the surf. pressure gradient trend
	91	!!
[508]	92	!! References : Madec et al. 1988, ocean modelling, issue 78, 1-6.
[3]	93	!!---------------------------------------------------------------------
	94	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[508]	95	INTEGER, INTENT( out ) :: kindic ! solver flag (<0 when the solver does not converge)
[3]	96	INTEGER :: ji, jj, jk ! dummy loop indices
[314]	97	REAL(wp) :: zbsfa, zgcx, z2dt
[3]	98	# if defined key_obc
	99	INTEGER :: ip, ii, ij
	100	INTEGER :: iii, ijj, jip, jnic
	101	INTEGER :: it, itm, itm2, ib, ibm, ibm2
	102	REAL(wp) :: z2dtr
	103	# endif
	104	!!----------------------------------------------------------------------
	105
	106	IF( kt == nit000 ) THEN
	107	IF(lwp) WRITE(numout,*)
	108	IF(lwp) WRITE(numout,*) 'dyn_spg_rl : surface pressure gradient trend'
	109	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
	110
	111	! set to zero rigid-lid specific arrays
[508]	112	spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction)
	113	spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction)
	114
	115	CALL solver_init( nit000 ) ! Elliptic solver initialisation
	116
	117	! read rigid lid arrays in restart file
	118	CALL rl_rst( nit000, 'READ' ) ! read or initialize the following fields:
	119	! ! gcx, gcxb, bsfb, bsfn, bsfd
[3]	120	ENDIF
	121
[508]	122	! Vertically averaged momentum trend
	123	! ------------------------------------
[3]	124	! ! ===============
	125	DO jj = 2, jpjm1 ! Vertical slab
	126	! ! ===============
[508]	127
	128	spgu(:,jj) = 0. ! initialization to zero
[3]	129	spgv(:,jj) = 0.
[508]	130	DO jk = 1, jpk ! vertical sum
[3]	131	DO ji = 2, jpim1
	132	spgu(ji,jj) = spgu(ji,jj) + ua(ji,jj,jk) * fse3u(ji,jj,jk) * umask(ji,jj,jk)
	133	spgv(ji,jj) = spgv(ji,jj) + va(ji,jj,jk) * fse3v(ji,jj,jk) * vmask(ji,jj,jk)
	134	END DO
	135	END DO
[508]	136	spgu(:,jj) = spgu(:,jj) * hur(:,jj) ! divide by the depth
[3]	137	spgv(:,jj) = spgv(:,jj) * hvr(:,jj)
	138
	139	! ! ===============
	140	END DO ! End of slab
	141	! ! ===============
	142
	143	! Boundary conditions on (spgu,spgv)
	144	CALL lbc_lnk( spgu, 'U', -1. )
	145	CALL lbc_lnk( spgv, 'V', -1. )
	146
[508]	147	! Barotropic streamfunction trend (bsfd)
[3]	148	! ----------------------------------
[508]	149	DO jj = 2, jpjm1 ! Curl of the vertically averaged velocity
[3]	150	DO ji = 2, jpim1
	151	gcb(ji,jj) = -gcdprc(ji,jj) &
[508]	152	& ( ( e2v(ji+1,jj )spgv(ji+1,jj ) - e2v(ji,jj)*spgv(ji,jj) ) &
	153	& -( e1u(ji ,jj+1)spgu(ji ,jj+1) - e1u(ji,jj)spgu(ji,jj) ) )
[3]	154	END DO
	155	END DO
	156
	157	# if defined key_obc
[508]	158	DO jj = 2, jpjm1 ! Open boundary contribution
[3]	159	DO ji = 2, jpim1
	160	gcb(ji,jj) = gcb(ji,jj) - gcdprc(ji,jj) * gcbob(ji,jj)
	161	END DO
	162	END DO
	163	# else
	164	! No open boundary contribution
	165	# endif
	166
	167	! First guess using previous solution of the elliptic system and
	168	! not bsfd since the system is solved with 0 as coastal boundary
	169	! condition. Also include a swap array (gcx,gxcb)
	170	DO jj = 2, jpjm1
	171	DO ji = 2, jpim1
	172	zgcx = gcx(ji,jj)
	173	gcx (ji,jj) = 2.*zgcx - gcxb(ji,jj)
	174	gcxb(ji,jj) = zgcx
	175	END DO
	176	END DO
[314]	177	! applied the lateral boundary conditions
[784]	178	IF( nsolv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 ) CALL lbc_lnk_e( gcb, c_solver_pt, 1. )
[3]	179
	180	! Relative precision (computation on one processor)
	181	rnorme = 0.e0
[314]	182	rnorme = SUM( gcb(1:nlci,1:nlcj) * gcdmat(1:nlci,1:nlcj) * gcb(1:nlci,1:nlcj) * bmask(:,:) )
[31]	183	IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain
	184
[3]	185	epsr = epsepsrnorme
	186	ncut = 0
	187	! if rnorme is 0, the solution is 0, the solver isn't called
	188	IF( rnorme == 0.e0 ) THEN
	189	bsfd (:,:) = 0.e0
	190	res = 0.e0
	191	niter = 0
	192	ncut = 999
	193	ENDIF
	194
	195	kindic = 0
	196	! solve the bsf system ===> solution in gcx array
	197	IF( ncut == 0 ) THEN
	198	SELECT CASE ( nsolv )
	199	CASE ( 1 ) ! diagonal preconditioned conjuguate gradient
	200	CALL sol_pcg( kindic )
	201	CASE( 2 ) ! successive-over-relaxation
[31]	202	CALL sol_sor( kindic )
[3]	203	CASE( 3 ) ! FETI solver
	204	CALL sol_fet( kindic )
	205	CASE DEFAULT ! e r r o r in nsolv namelist parameter
[474]	206	WRITE(ctmp1,*) ' ~~~~~~~~~~ not = ', nsolv
[784]	207	CALL ctl_stop( ' dyn_spg_rl : e r r o r, nsolv = 1, 2 or 3', ctmp1 )
[3]	208	END SELECT
	209	ENDIF
	210
	211	! bsf trend update
	212	! ----------------
[314]	213	bsfd(1:nlci,1:nlcj) = gcx(1:nlci,1:nlcj)
[3]	214
	215	! update bsf trend with islands trend
	216	! -----------------------------------
[78]	217	IF( lk_isl ) CALL isl_dyn_spg ! update bsfd
[3]	218
	219	# if defined key_obc
	220	! Compute bsf trend for OBC points (not masked)
	221
[78]	222	IF( lp_obc_east ) THEN
[3]	223	! compute bsf trend at the boundary from bsfeob, computed in obc_spg
	224	IF( neuler == 0 .AND. kt == nit000 ) THEN
	225	z2dtr = 1. / rdt
	226	ELSE
	227	z2dtr = 1. / (2. * rdt )
	228	ENDIF
	229	! (jped,jpefm1),nieob
	230	DO ji = fs_nie0, fs_nie1 ! vector opt.
	231	DO jj = nje0m1, nje1m1
	232	bsfd(ji,jj) = ( bsfeob(jj) - bsfb(ji,jj) ) * z2dtr
	233	END DO
	234	END DO
	235	ENDIF
	236
[78]	237	IF( lp_obc_west ) THEN
[3]	238	! compute bsf trend at the boundary from bsfwob, computed in obc_spg
	239	IF( neuler == 0 .AND. kt == nit000 ) THEN
	240	z2dtr = 1. / rdt
	241	ELSE
	242	z2dtr = 1. / ( 2. * rdt )
	243	ENDIF
	244	! (jpwd,jpwfm1),niwob
	245	DO ji = fs_niw0, fs_niw1 ! vector opt.
	246	DO jj = njw0m1, njw1m1
	247	bsfd(ji,jj) = ( bsfwob(jj) - bsfb(ji,jj) ) * z2dtr
	248	END DO
	249	END DO
	250	ENDIF
	251
[78]	252	IF( lp_obc_north ) THEN
[3]	253	! compute bsf trend at the boundary from bsfnob, computed in obc_spg
	254	IF( neuler == 0 .AND. kt == nit000 ) THEN
	255	z2dtr = 1. / rdt
	256	ELSE
	257	z2dtr = 1. / ( 2. * rdt )
	258	ENDIF
	259	! njnob,(jpnd,jpnfm1)
	260	DO jj = fs_njn0, fs_njn1 ! vector opt.
	261	DO ji = nin0m1, nin1m1
	262	bsfd(ji,jj) = ( bsfnob(ji) - bsfb(ji,jj) ) * z2dtr
	263	END DO
	264	END DO
	265	ENDIF
	266
[78]	267	IF( lp_obc_south ) THEN
[3]	268	! compute bsf trend at the boundary from bsfsob, computed in obc_spg
	269	IF( neuler == 0 .AND. kt == nit000 ) THEN
	270	z2dtr = 1. / rdt
	271	ELSE
	272	z2dtr = 1. / ( 2. * rdt )
	273	ENDIF
	274	! njsob,(jpsd,jpsfm1)
	275	DO jj = fs_njs0, fs_njs1 ! vector opt.
	276	DO ji = nis0m1, nis1m1
	277	bsfd(ji,jj) = ( bsfsob(ji) - bsfb(ji,jj) ) * z2dtr
	278	END DO
	279	END DO
	280	ENDIF
	281
	282	! compute bsf trend for isolated coastline points
	283
	284	IF( neuler == 0 .AND. kt == nit000 ) THEN
	285	z2dtr = 1. / rdt
	286	ELSE
	287	z2dtr = 1. /( 2. * rdt )
	288	ENDIF
	289
	290	IF( nbobc > 1 ) THEN
	291	DO jnic = 1,nbobc - 1
	292	ip = mnic(0,jnic)
	293	DO jip = 1,ip
	294	ii = miic(jip,0,jnic)
	295	ij = mjic(jip,0,jnic)
	296	IF( ii >= 1 + nimpp - 1 .AND. ii <= jpi + nimpp -1 .AND. &
	297	ij >= 1 + njmpp - 1 .AND. ij <= jpj + njmpp -1 ) THEN
	298	iii = ii - nimpp + 1
	299	ijj = ij - njmpp + 1
	300	bsfd(iii,ijj) = ( bsfic(jnic) - bsfb(iii,ijj) ) * z2dtr
	301	ENDIF
	302	END DO
	303	END DO
	304	ENDIF
	305	# endif
	306
	307	! 4. Barotropic stream function and array swap
	308	! --------------------------------------------
	309	! Leap-frog time scheme, time filter and array swap
	310	IF( neuler == 0 .AND. kt == nit000 ) THEN
	311	! Euler time stepping (first time step, starting from rest)
	312	z2dt = rdt
	313	DO jj = 1, jpj
	314	DO ji = 1, jpi
	315	zbsfa = bsfb(ji,jj) + z2dt * bsfd(ji,jj)
	316	bsfb(ji,jj) = bsfn(ji,jj)
	317	bsfn(ji,jj) = zbsfa
	318	END DO
	319	END DO
	320	ELSE
	321	! Leap-frog time stepping - Asselin filter
	322	z2dt = 2.*rdt
	323	DO jj = 1, jpj
	324	DO ji = 1, jpi
	325	zbsfa = bsfb(ji,jj) + z2dt * bsfd(ji,jj)
	326	bsfb(ji,jj) = atfp * ( bsfb(ji,jj) + zbsfa ) + atfp1 * bsfn(ji,jj)
	327	bsfn(ji,jj) = zbsfa
	328	END DO
	329	END DO
	330	ENDIF
	331
	332	# if defined key_obc
[508]	333	ib = 1 ! space index on boundary arrays
	334	ibm = 2
	335	ibm2 = 3
	336	it = 1 ! time index on boundary arrays
	337	itm = 2
	338	itm2 = 3
	339
[3]	340	! Swap of boundary arrays
[78]	341	IF( lp_obc_east ) THEN
[3]	342	! (jped,jpef),nieob
	343	IF( kt < nit000+3 .AND. .NOT.ln_rstart ) THEN
	344	DO jj = nje0m1, nje1
	345	! fields itm2 <== itm
	346	bebnd(jj,ib ,itm2) = bebnd(jj,ib ,itm)
	347	bebnd(jj,ibm ,itm2) = bebnd(jj,ibm ,itm)
	348	bebnd(jj,ibm2,itm2) = bebnd(jj,ibm2,itm)
	349	bebnd(jj,ib ,itm ) = bebnd(jj,ib ,it )
	350	END DO
	351	ELSE
	352	! fields itm <== it plus time filter at the boundary
	353	DO ji = fs_nie0, fs_nie1 ! vector opt.
	354	DO jj = nje0m1, nje1
	355	bebnd(jj,ib ,itm2) = bebnd(jj,ib ,itm)
	356	bebnd(jj,ibm ,itm2) = bebnd(jj,ibm ,itm)
	357	bebnd(jj,ibm2,itm2) = bebnd(jj,ibm2,itm)
	358	bebnd(jj,ib ,itm ) = atfp * ( bebnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bebnd(jj,ib,it)
	359	bebnd(jj,ibm ,itm ) = bebnd(jj,ibm ,it )
	360	bebnd(jj,ibm2,itm ) = bebnd(jj,ibm2,it )
	361	END DO
	362	END DO
	363	ENDIF
	364	! fields it <== now (kt+1)
	365	DO ji = fs_nie0, fs_nie1 ! vector opt.
	366	DO jj = nje0m1, nje1
	367	bebnd(jj,ib ,it ) = bsfn (ji ,jj)
	368	bebnd(jj,ibm ,it ) = bsfn (ji-1,jj)
	369	bebnd(jj,ibm2,it ) = bsfn (ji-2,jj)
	370	END DO
	371	END DO
[31]	372	IF( lk_mpp ) CALL mppobc( bebnd, jpjed, jpjef, jpieob, 3*3, 2, jpj )
[3]	373	ENDIF
	374
[78]	375	IF( lp_obc_west ) THEN
[3]	376	! (jpwd,jpwf),niwob
	377	IF( kt < nit000+3 .AND. .NOT.ln_rstart ) THEN
	378	DO jj = njw0m1, njw1
	379	! fields itm2 <== itm
	380	bwbnd(jj,ib ,itm2) = bwbnd(jj,ib ,itm)
	381	bwbnd(jj,ibm ,itm2) = bwbnd(jj,ibm ,itm)
	382	bwbnd(jj,ibm2,itm2) = bwbnd(jj,ibm2,itm)
	383	bwbnd(jj,ib ,itm ) = bwbnd(jj,ib ,it )
	384	END DO
	385	ELSE
	386	DO ji = fs_niw0, fs_niw1 ! Vector opt.
	387	DO jj = njw0m1, njw1
	388	bwbnd(jj,ib ,itm2) = bwbnd(jj,ib ,itm)
	389	bwbnd(jj,ibm ,itm2) = bwbnd(jj,ibm ,itm)
	390	bwbnd(jj,ibm2,itm2) = bwbnd(jj,ibm2,itm)
	391	! fields itm <== it plus time filter at the boundary
	392	bwbnd(jj,ib ,itm ) = atfp * ( bwbnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bwbnd(jj,ib,it)
	393	bwbnd(jj,ibm ,itm ) = bwbnd(jj,ibm ,it )
	394	bwbnd(jj,ibm2,itm ) = bwbnd(jj,ibm2,it )
	395	END DO
	396	END DO
	397	ENDIF
	398	! fields it <== now (kt+1)
	399	DO ji = fs_niw0, fs_niw1 ! Vector opt.
	400	DO jj = njw0m1, njw1
	401	bwbnd(jj,ib ,it ) = bsfn (ji ,jj)
	402	bwbnd(jj,ibm ,it ) = bsfn (ji+1,jj)
	403	bwbnd(jj,ibm2,it ) = bsfn (ji+2,jj)
	404	END DO
	405	END DO
[31]	406	IF( lk_mpp ) CALL mppobc( bwbnd, jpjwd, jpjwf, jpiwob, 3*3, 2, jpj )
[3]	407	ENDIF
	408
[78]	409	IF( lp_obc_north ) THEN
[3]	410	! njnob,(jpnd,jpnf)
	411	IF( kt < nit000 + 3 .AND. .NOT.ln_rstart ) THEN
	412	DO ji = nin0m1, nin1
	413	! fields itm2 <== itm
	414	bnbnd(ji,ib ,itm2) = bnbnd(ji,ib ,itm)
	415	bnbnd(ji,ibm ,itm2) = bnbnd(ji,ibm ,itm)
	416	bnbnd(ji,ibm2,itm2) = bnbnd(ji,ibm2,itm)
	417	bnbnd(ji,ib ,itm ) = bnbnd(ji,ib ,it )
	418	END DO
	419	ELSE
	420	DO jj = fs_njn0, fs_njn1 ! Vector opt.
	421	DO ji = nin0m1, nin1
	422	bnbnd(ji,ib ,itm2) = bnbnd(ji,ib ,itm)
	423	bnbnd(ji,ibm ,itm2) = bnbnd(ji,ibm ,itm)
	424	bnbnd(ji,ibm2,itm2) = bnbnd(ji,ibm2,itm)
	425	! fields itm <== it plus time filter at the boundary
	426	bnbnd(jj,ib ,itm ) = atfp * ( bnbnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bnbnd(jj,ib,it)
	427	bnbnd(ji,ibm ,itm ) = bnbnd(ji,ibm ,it )
	428	bnbnd(ji,ibm2,itm ) = bnbnd(ji,ibm2,it )
	429	END DO
	430	END DO
	431	ENDIF
	432	! fields it <== now (kt+1)
	433	DO jj = fs_njn0, fs_njn1 ! Vector opt.
	434	DO ji = nin0m1, nin1
	435	bnbnd(ji,ib ,it ) = bsfn (ji,jj )
	436	bnbnd(ji,ibm ,it ) = bsfn (ji,jj-1)
	437	bnbnd(ji,ibm2,it ) = bsfn (ji,jj-2)
	438	END DO
	439	END DO
[31]	440	IF( lk_mpp ) CALL mppobc( bnbnd, jpind, jpinf, jpjnob, 3*3, 1, jpi )
[3]	441	ENDIF
	442
[78]	443	IF( lp_obc_south ) THEN
[31]	444	! njsob,(jpsd,jpsf)
	445	IF( kt < nit000+3 .AND. .NOT.ln_rstart ) THEN
[3]	446	DO ji = nis0m1, nis1
[31]	447	! fields itm2 <== itm
[3]	448	bsbnd(ji,ib ,itm2) = bsbnd(ji,ib ,itm)
	449	bsbnd(ji,ibm ,itm2) = bsbnd(ji,ibm ,itm)
	450	bsbnd(ji,ibm2,itm2) = bsbnd(ji,ibm2,itm)
[31]	451	bsbnd(ji,ib ,itm ) = bsbnd(ji,ib ,it )
[3]	452	END DO
[31]	453	ELSE
	454	DO jj = fs_njs0, fs_njs1 ! vector opt.
	455	DO ji = nis0m1, nis1
	456	bsbnd(ji,ib ,itm2) = bsbnd(ji,ib ,itm)
	457	bsbnd(ji,ibm ,itm2) = bsbnd(ji,ibm ,itm)
	458	bsbnd(ji,ibm2,itm2) = bsbnd(ji,ibm2,itm)
	459	! fields itm <== it plus time filter at the boundary
	460	bsbnd(jj,ib ,itm ) = atfp * ( bsbnd(jj,ib,itm) + bsfn(ji,jj) ) + atfp1 * bsbnd(jj,ib,it)
	461	bsbnd(ji,ibm ,itm ) = bsbnd(ji,ibm ,it )
	462	bsbnd(ji,ibm2,itm ) = bsbnd(ji,ibm2,it )
	463	END DO
	464	END DO
	465	ENDIF
	466	DO jj = fs_njs0, fs_njs1 ! vector opt.
	467	DO ji = nis0m1, nis1
	468	! fields it <== now (kt+1)
	469	bsbnd(ji,ib ,it ) = bsfn (ji,jj )
	470	bsbnd(ji,ibm ,it ) = bsfn (ji,jj+1)
	471	bsbnd(ji,ibm2,it ) = bsfn (ji,jj+2)
	472	END DO
[3]	473	END DO
[31]	474	IF( lk_mpp ) CALL mppobc( bsbnd, jpisd, jpisf, jpjsob, 3*3, 1, jpi )
[3]	475	ENDIF
	476	# endif
	477
	478	! add the surface pressure trend to the general trend
	479	! -----------------------------------------------------
	480	DO jj = 2, jpjm1
	481	! update the surface pressure gradient with the barotropic trend
	482	DO ji = 2, jpim1
	483	spgu(ji,jj) = spgu(ji,jj) + hur(ji,jj) / e2u(ji,jj) * ( bsfd(ji,jj) - bsfd(ji ,jj-1) )
	484	spgv(ji,jj) = spgv(ji,jj) - hvr(ji,jj) / e1v(ji,jj) * ( bsfd(ji,jj) - bsfd(ji-1,jj ) )
	485	END DO
	486	! add the surface pressure gradient trend to the general trend
	487	DO jk = 1, jpkm1
	488	DO ji = 2, jpim1
	489	ua(ji,jj,jk) = ua(ji,jj,jk) - spgu(ji,jj)
	490	va(ji,jj,jk) = va(ji,jj,jk) - spgv(ji,jj)
	491	END DO
	492	END DO
	493	END DO
	494
[508]	495	! write rigid lid arrays in restart file
	496	! --------------------------------------
	497	IF( lrst_oce ) CALL rl_rst( kt, 'WRITE' )
	498
[3]	499	END SUBROUTINE dyn_spg_rl
	500
[508]	501
	502	SUBROUTINE rl_rst( kt, cdrw )
	503	!!---------------------------------------------------------------------
	504	!! * ROUTINE rl_rst *
	505	!!
	506	!! ** Purpose : Read or write rigid-lid arrays in ocean restart file
	507	!!----------------------------------------------------------------------
	508	INTEGER , INTENT(in) :: kt ! ocean time-step
	509	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	510	!!----------------------------------------------------------------------
	511	!
	512	IF( TRIM(cdrw) == 'READ' ) THEN
[746]	513	IF( iom_varid( numror, 'gcx', ldstop = .FALSE. ) > 0 ) THEN
[508]	514	! Caution : extra-hallow
	515	! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
[683]	516	CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
	517	CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
	518	CALL iom_get( numror, jpdom_autoglo, 'bsfb', bsfb(:,:) )
	519	CALL iom_get( numror, jpdom_autoglo, 'bsfn', bsfn(:,:) )
	520	CALL iom_get( numror, jpdom_autoglo, 'bsfd', bsfd(:,:) )
[508]	521	IF( neuler == 0 ) THEN
	522	gcxb(:,:) = gcx (:,:)
	523	bsfb(:,:) = bsfn(:,:)
	524	ENDIF
	525	ELSE
	526	gcx (:,:) = 0.e0
	527	gcxb(:,:) = 0.e0
	528	bsfb(:,:) = 0.e0
	529	bsfn(:,:) = 0.e0
	530	bsfd(:,:) = 0.e0
	531	ENDIF
	532	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
	533	! Caution : extra-hallow, gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
	534	CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
	535	CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
	536	CALL iom_rstput( kt, nitrst, numrow, 'bsfb', bsfb(:,:) )
	537	CALL iom_rstput( kt, nitrst, numrow, 'bsfn', bsfn(:,:) )
	538	CALL iom_rstput( kt, nitrst, numrow, 'bsfd', bsfd(:,:) )
	539	ENDIF
	540	!
	541	END SUBROUTINE rl_rst
	542
[3]	543	#else
	544	!!----------------------------------------------------------------------
	545	!! 'key_dynspg_rl' NO rigid lid
	546	!!----------------------------------------------------------------------
	547	CONTAINS
	548	SUBROUTINE dyn_spg_rl( kt, kindic ) ! Empty routine
[31]	549	WRITE(,) 'dyn_spg_rl: You should not have seen this print! error?', kt, kindic
[3]	550	END SUBROUTINE dyn_spg_rl
	551	#endif
	552
	553	!!======================================================================
	554	END MODULE dynspg_rl

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

Download in other formats: