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.F90 in branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90 @ 3957

Last change on this file since 3957 was 3957, checked in by acc, 11 years ago
Branch 2013/dev_r3858_NOC_ZTC, #863. Fixes to tidal potential forcing changes to ensure no numerical changes due to the re-styling
Property svn:keywords set to `Id`
File size: 13.6 KB

Line
1	MODULE dynspg
2	!!======================================================================
3	!! * MODULE dynspg *
4	!! Ocean dynamics: surface pressure gradient control
5	!!======================================================================
6	!! History : 1.0 ! 2005-12 (C. Talandier, G. Madec, V. Garnier) Original code
7	!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
8	!!----------------------------------------------------------------------
9
10	!!----------------------------------------------------------------------
11	!! dyn_spg : update the dynamics trend with the lateral diffusion
12	!! dyn_spg_ctl : initialization, namelist read, and parameters control
13	!!----------------------------------------------------------------------
14	USE oce ! ocean dynamics and tracers variables
15	USE dom_oce ! ocean space and time domain variables
16	USE phycst ! physical constants
17	USE sbc_oce ! surface boundary condition: ocean
18	USE sbcapr ! surface boundary condition: atmospheric pressure
19	USE dynspg_oce ! surface pressure gradient variables
20	USE dynspg_exp ! surface pressure gradient (dyn_spg_exp routine)
21	USE dynspg_ts ! surface pressure gradient (dyn_spg_ts routine)
22	USE dynspg_flt ! surface pressure gradient (dyn_spg_flt routine)
23	USE dynadv ! dynamics: vector invariant versus flux form
24	USE sbctide
25	USE updtide
26	USE trdmod ! ocean dynamics trends
27	USE trdmod_oce ! ocean variables trends
28	USE prtctl ! Print control (prt_ctl routine)
29	USE in_out_manager ! I/O manager
30	USE lib_mpp ! MPP library
31	USE solver ! solver initialization
32	USE wrk_nemo ! Memory Allocation
33	USE timing ! Timing
34
35
36	IMPLICIT NONE
37	PRIVATE
38
39	PUBLIC dyn_spg ! routine called by step module
40	PUBLIC dyn_spg_init ! routine called by opa module
41
42	INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from lk_dynspg_...
43
44	!! * Substitutions
45	# include "domzgr_substitute.h90"
46	# include "vectopt_loop_substitute.h90"
47	!!----------------------------------------------------------------------
48	!! NEMO/OPA 3.2 , LODYC-IPSL (2009)
49	!! $Id$
50	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
51	!!----------------------------------------------------------------------
52	CONTAINS
53
54	SUBROUTINE dyn_spg( kt, kindic )
55	!!----------------------------------------------------------------------
56	!! * ROUTINE dyn_spg *
57	!!
58	!! ** Purpose : achieve the momentum time stepping by computing the
59	!! last trend, the surface pressure gradient including the
60	!! atmospheric pressure forcing (ln_apr_dyn=T), and performing
61	!! the Leap-Frog integration.
62	!!gm In the current version only the filtered solution provide
63	!!gm the after velocity, in the 2 other (ua,va) are still the trends
64	!!
65	!! ** Method : Three schemes:
66	!! - explicit computation : the spg is evaluated at now
67	!! - filtered computation : the Roulet & madec (2000) technique is used
68	!! - split-explicit computation: a time splitting technique is used
69	!!
70	!! ln_apr_dyn=T : the atmospheric pressure forcing is applied
71	!! as the gradient of the inverse barometer ssh:
72	!! apgu = - 1/rau0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb]
73	!! apgv = - 1/rau0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb]
74	!! Note that as all external forcing a time averaging over a two rdt
75	!! period is used to prevent the divergence of odd and even time step.
76	!!
77	!! N.B. : When key_esopa is used all the scheme are tested, regardless
78	!! of the physical meaning of the results.
79	!!----------------------------------------------------------------------
80	!
81	INTEGER, INTENT(in ) :: kt ! ocean time-step index
82	INTEGER, INTENT( out) :: kindic ! solver flag
83	!
84	INTEGER :: ji, jj, jk ! dummy loop indices
85	REAL(wp) :: z2dt, zg_2, zintp, zgrau0r ! temporary scalar
86	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
87	REAL(wp), POINTER, DIMENSION(:,:) :: zpice
88	!!----------------------------------------------------------------------
89	!
90	IF( nn_timing == 1 ) CALL timing_start('dyn_spg')
91	!
92
93	!!gm NOTA BENE : the dynspg_exp and dynspg_ts should be modified so that
94	!!gm they return the after velocity, not the trends (as in trazdf_imp...)
95	!!gm In this case, change/simplify dynnxt
96
97
98	IF( l_trddyn ) THEN ! temporary save of ta and sa trends
99	CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )
100	ztrdu(:,:,:) = ua(:,:,:)
101	ztrdv(:,:,:) = va(:,:,:)
102	ENDIF
103
104	IF( ln_apr_dyn & ! atmos. pressure
105	.OR. ( .NOT.lk_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) ) & ! tide potential (no time slitting)
106	.OR. nn_ice_embd == 2 ) THEN ! embedded sea-ice
107	!
108	DO jj = 2, jpjm1
109	DO ji = fs_2, fs_jpim1 ! vector opt.
110	spgu(ji,jj) = 0._wp
111	spgv(ji,jj) = 0._wp
112	END DO
113	END DO
114	!
115	IF( ln_apr_dyn ) THEN !== Atmospheric pressure gradient ==!
116	zg_2 = grav * 0.5
117	DO jj = 2, jpjm1 ! gradient of Patm using inverse barometer ssh
118	DO ji = fs_2, fs_jpim1 ! vector opt.
119	spgu(ji,jj) = spgu(ji,jj) + zg_2 * ( ssh_ib (ji+1,jj) - ssh_ib (ji,jj) &
120	& + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj) ) /e1u(ji,jj)
121	spgv(ji,jj) = spgv(ji,jj) + zg_2 * ( ssh_ib (ji,jj+1) - ssh_ib (ji,jj) &
122	& + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)
123	END DO
124	END DO
125	ENDIF
126	!
127	! !== tide potential forcing term ==!
128	IF( .NOT.lk_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case
129	!
130	CALL upd_tide( kt ) ! update tide potential
131	!
132	DO jj = 2, jpjm1 ! add tide potential forcing
133	DO ji = fs_2, fs_jpim1 ! vector opt.
134	spgv(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
135	spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
136	END DO
137	END DO
138	ENDIF
139	!
140	IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
141	CALL wrk_alloc( jpi, jpj, zpice )
142	!
143	zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
144	zgrau0r = - grav * r1_rau0
145	zpice(:,:) = ( zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:) ) * zgrau0r
146	DO jj = 2, jpjm1
147	DO ji = fs_2, fs_jpim1 ! vector opt.
148	spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
149	spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
150	END DO
151	END DO
152	!
153	CALL wrk_dealloc( jpi, jpj, zpice )
154	ENDIF
155	!
156	DO jk = 1, jpkm1 !== Add all terms to the general trend
157	DO jj = 2, jpjm1
158	DO ji = fs_2, fs_jpim1 ! vector opt.
159	ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
160	va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
161	END DO
162	END DO
163	END DO
164	ENDIF
165
166	SELECT CASE ( nspg ) ! compute surf. pressure gradient trend and add it to the general trend
167	!
168	CASE ( 0 ) ; CALL dyn_spg_exp( kt ) ! explicit
169	CASE ( 1 ) ; CALL dyn_spg_ts ( kt ) ! time-splitting
170	CASE ( 2 ) ; CALL dyn_spg_flt( kt, kindic ) ! filtered
171	!
172	CASE ( -1 ) ! esopa: test all possibility with control print
173	CALL dyn_spg_exp( kt )
174	CALL prt_ctl( tab3d_1=ua, clinfo1=' spg0 - Ua: ', mask1=umask, &
175	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
176	CALL dyn_spg_ts ( kt )
177	CALL prt_ctl( tab3d_1=ua, clinfo1=' spg1 - Ua: ', mask1=umask, &
178	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
179	CALL dyn_spg_flt( kt, kindic )
180	CALL prt_ctl( tab3d_1=ua, clinfo1=' spg2 - Ua: ', mask1=umask, &
181	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
182	END SELECT
183	!
184	IF( l_trddyn ) THEN ! save the surface pressure gradient trends for further diagnostics
185	SELECT CASE ( nspg )
186	CASE ( 0, 1 )
187	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
188	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
189	CASE( 2 )
190	z2dt = 2. * rdt
191	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
192	ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:)
193	ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:)
194	END SELECT
195	CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_spg, 'DYN', kt )
196	!
197	CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )
198	ENDIF
199	! ! print mean trends (used for debugging)
200	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' spg - Ua: ', mask1=umask, &
201	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
202	!
203	IF( nn_timing == 1 ) CALL timing_stop('dyn_spg')
204	!
205	END SUBROUTINE dyn_spg
206
207
208	SUBROUTINE dyn_spg_init
209	!!---------------------------------------------------------------------
210	!! * ROUTINE dyn_spg_init *
211	!!
212	!! ** Purpose : Control the consistency between cpp options for
213	!! surface pressure gradient schemes
214	!!----------------------------------------------------------------------
215	INTEGER :: ioptio
216	!!----------------------------------------------------------------------
217	!
218	IF( nn_timing == 1 ) CALL timing_start('dyn_spg_init')
219	!
220	IF(lwp) THEN ! Control print
221	WRITE(numout,*)
222	WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme'
223	WRITE(numout,*) '~~~~~~~~~~~'
224	WRITE(numout,*) ' Explicit free surface lk_dynspg_exp = ', lk_dynspg_exp
225	WRITE(numout,*) ' Free surface with time splitting lk_dynspg_ts = ', lk_dynspg_ts
226	WRITE(numout,*) ' Filtered free surface cst volume lk_dynspg_flt = ', lk_dynspg_flt
227	ENDIF
228
229	! ! allocate dyn_spg arrays
230	IF( lk_dynspg_ts ) THEN
231	IF( dynspg_oce_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_oce arrays')
232	IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays')
233	ENDIF
234
235	! ! Control of surface pressure gradient scheme options
236	ioptio = 0
237	IF(lk_dynspg_exp) ioptio = ioptio + 1
238	IF(lk_dynspg_ts ) ioptio = ioptio + 1
239	IF(lk_dynspg_flt) ioptio = ioptio + 1
240	!
241	IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ioptio == 0 ) &
242	& CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' )
243	!
244	IF( lk_esopa ) nspg = -1
245	IF( lk_dynspg_exp) nspg = 0
246	IF( lk_dynspg_ts ) nspg = 1
247	IF( lk_dynspg_flt) nspg = 2
248	!
249	IF( lk_esopa ) nspg = -1
250	!
251	IF(lwp) THEN
252	WRITE(numout,*)
253	IF( nspg == -1 ) WRITE(numout,*) ' ESOPA test All scheme used'
254	IF( nspg == 0 ) WRITE(numout,*) ' explicit free surface'
255	IF( nspg == 1 ) WRITE(numout,*) ' free surface with time splitting scheme'
256	IF( nspg == 2 ) WRITE(numout,*) ' filtered free surface'
257	ENDIF
258
259	#if defined key_dynspg_flt \|\| defined key_esopa
260	CALL solver_init( nit000 ) ! Elliptic solver initialisation
261	#endif
262
263	! ! Control of timestep choice
264	IF( lk_dynspg_ts .OR. lk_dynspg_exp ) THEN
265	IF( nn_cla == 1 ) CALL ctl_stop( 'Crossland advection not implemented for this free surface formulation' )
266	ENDIF
267
268	! ! Control of momentum formulation
269	IF( lk_dynspg_ts .AND. lk_vvl ) THEN
270	IF( .NOT.ln_dynadv_vec ) CALL ctl_stop( 'Flux form not implemented for this free surface formulation' )
271	ENDIF
272	!
273	IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_init')
274	!
275	END SUBROUTINE dyn_spg_init
276
277	!!======================================================================
278	END MODULE dynspg

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

Download in other formats: