Context Navigation

Back to Ticket #1126

Ticket #1126: istate.F90

File istate.F90, 25.5 KB (added by poddo, 11 years ago)

Line
1	MODULE istate
2	!!======================================================================
3	!! * MODULE istate *
4	!! Ocean state : initial state setting
5	!!=====================================================================
6	!! History : OPA ! 1989-12 (P. Andrich) Original code
7	!! 5.0 ! 1991-11 (G. Madec) rewritting
8	!! 6.0 ! 1996-01 (G. Madec) terrain following coordinates
9	!! 8.0 ! 2001-09 (M. Levy, M. Ben Jelloul) istate_eel
10	!! 8.0 ! 2001-09 (M. Levy, M. Ben Jelloul) istate_uvg
11	!! NEMO 1.0 ! 2003-08 (G. Madec, C. Talandier) F90: Free form, modules + EEL R5
12	!! - ! 2004-05 (A. Koch-Larrouy) istate_gyre
13	!! 2.0 ! 2006-07 (S. Masson) distributed restart using iom
14	!! 3.3 ! 2010-10 (C. Ethe) merge TRC-TRA
15	!! 3.4 ! 2011-04 (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn
16	!!----------------------------------------------------------------------
17
18	!!----------------------------------------------------------------------
19	!! istate_init : initial state setting
20	!! istate_tem : analytical profile for initial Temperature
21	!! istate_sal : analytical profile for initial Salinity
22	!! istate_eel : initial state setting of EEL R5 configuration
23	!! istate_gyre : initial state setting of GYRE configuration
24	!! istate_uvg : initial velocity in geostropic balance
25	!!----------------------------------------------------------------------
26	USE oce ! ocean dynamics and active tracers
27	USE dom_oce ! ocean space and time domain
28	USE daymod ! calendar
29	USE eosbn2 ! eq. of state, Brunt Vaisala frequency (eos routine)
30	USE ldftra_oce ! ocean active tracers: lateral physics
31	USE zdf_oce ! ocean vertical physics
32	USE phycst ! physical constants
33	USE dtatsd ! data temperature and salinity (dta_tsd routine)
34	USE restart ! ocean restart (rst_read routine)
35	USE in_out_manager ! I/O manager
36	USE iom ! I/O library
37	USE zpshde ! partial step: hor. derivative (zps_hde routine)
38	USE eosbn2 ! equation of state (eos bn2 routine)
39	USE domvvl ! varying vertical mesh
40	USE dynspg_oce ! pressure gradient schemes
41	USE dynspg_flt ! pressure gradient schemes
42	USE dynspg_exp ! pressure gradient schemes
43	USE dynspg_ts ! pressure gradient schemes
44	USE sol_oce ! ocean solver variables
45	USE lib_mpp ! MPP library
46	USE wrk_nemo ! Memory allocation
47	USE timing ! Timing
48
49	IMPLICIT NONE
50	PRIVATE
51
52	PUBLIC istate_init ! routine called by step.F90
53
54	!! * Substitutions
55	# include "domzgr_substitute.h90"
56	# include "vectopt_loop_substitute.h90"
57	!!----------------------------------------------------------------------
58	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
59	!! $Id: istate.F90 3711 2012-11-28 12:02:09Z flavoni $
60	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
61	!!----------------------------------------------------------------------
62	CONTAINS
63
64	SUBROUTINE istate_init
65	!!----------------------------------------------------------------------
66	!! * ROUTINE istate_init *
67	!!
68	!! ** Purpose : Initialization of the dynamics and tracer fields.
69	!!----------------------------------------------------------------------
70	! - ML - needed for initialization of e3t_b
71	INTEGER :: jk ! dummy loop indice
72	!!----------------------------------------------------------------------
73	!
74	IF( nn_timing == 1 ) CALL timing_start('istate_init')
75	!
76
77	IF(lwp) WRITE(numout,*)
78	IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
79	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
80
81	CALL dta_tsd_init ! Initialisation of T & S input data
82
83	rhd (:,:,: ) = 0.e0
84	rhop (:,:,: ) = 0.e0
85	rn2 (:,:,: ) = 0.e0
86	tsa (:,:,:,:) = 0.e0
87
88	IF( ln_rstart ) THEN ! Restart from a file
89	! ! -------------------
90	neuler = 1 ! Set time-step indicator at nit000 (leap-frog)
91	CALL rst_read ! Read the restart file
92	! ! define e3u_b, e3v_b from e3t_b read in restart file
93	CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) )
94	CALL day_init ! model calendar (using both namelist and restart infos)
95	ELSE
96	! ! Start from rest
97	! ! ---------------
98	numror = 0 ! define numror = 0 -> no restart file to read
99	neuler = 0 ! Set time-step indicator at nit000 (euler forward)
100	CALL day_init ! model calendar (using both namelist and restart infos)
101	! ! Initialization of ocean to zero
102	! before fields ! now fields
103	sshb (:,:) = 0._wp ; sshn (:,:) = 0._wp
104	ub (:,:,:) = 0._wp ; un (:,:,:) = 0._wp
105	vb (:,:,:) = 0._wp ; vn (:,:,:) = 0._wp
106	rotb (:,:,:) = 0._wp ; rotn (:,:,:) = 0._wp
107	hdivb(:,:,:) = 0._wp ; hdivn(:,:,:) = 0._wp
108	!
109	IF( cp_cfg == 'eel' ) THEN
110	CALL istate_eel ! EEL configuration : start from pre-defined U,V T-S fields
111	ELSEIF( cp_cfg == 'gyre' ) THEN
112	CALL istate_gyre ! GYRE configuration : start from pre-defined T-S fields
113	ELSEIF( ln_tsd_init ) THEN ! Initial T-S fields read in files
114	CALL dta_tsd( nit000, tsb ) ! read 3D T and S data at nit000
115	tsn(:,:,:,:) = tsb(:,:,:,:)
116	!
117	ELSE ! Initial T-S fields defined analytically
118	CALL istate_t_s
119	ENDIF
120	!
121	CALL eos( tsb, rhd, rhop ) ! before potential and in situ densities
122	#if ! defined key_c1d
123	IF( ln_zps ) CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv, & ! zps: before hor. gradient
124	& rhd, gru , grv ) ! of t,s,rd at ocean bottom
125	#endif
126	!
127	! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
128	IF( lk_vvl ) THEN
129	DO jk = 1, jpk
130	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
131	ENDDO
132	ENDIF
133	! ! define e3u_b, e3v_b from e3t_b initialized in domzgr
134	CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) )
135	!
136	ENDIF
137	!
138	IF( lk_agrif ) THEN ! read free surface arrays in restart file
139	IF( ln_rstart ) THEN
140	IF( lk_dynspg_flt ) THEN ! read or initialize the following fields
141	! ! gcx, gcxb for agrif_opa_init
142	IF( sol_oce_alloc() > 0 ) CALL ctl_stop('agrif sol_oce_alloc: allocation of arrays failed')
143	CALL flt_rst( nit000, 'READ' )
144	ENDIF
145	ENDIF ! explicit case not coded yet with AGRIF
146	ENDIF
147	!
148	IF( nn_timing == 1 ) CALL timing_stop('istate_init')
149	!
150	END SUBROUTINE istate_init
151
152	SUBROUTINE istate_t_s
153	!!---------------------------------------------------------------------
154	!! * ROUTINE istate_t_s *
155	!!
156	!! ** Purpose : Intialization of the temperature field with an
157	!! analytical profile or a file (i.e. in EEL configuration)
158	!!
159	!! ** Method : - temperature: use Philander analytic profile
160	!! - salinity : use to a constant value 35.5
161	!!
162	!! References : Philander ???
163	!!----------------------------------------------------------------------
164	INTEGER :: ji, jj, jk
165	REAL(wp) :: zsal = 35.50
166	!!----------------------------------------------------------------------
167	!
168	IF(lwp) WRITE(numout,*)
169	IF(lwp) WRITE(numout,*) 'istate_t_s : Philander s initial temperature profile'
170	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ and constant salinity (',zsal,' psu)'
171	!
172	! DO jk = 1, jpk
173	! tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) ) &
174	! & + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.) ) * tmask(:,:,jk)
175	! tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
176	! END DO
177	tsn(:,:,:,jp_tem) = 15.0 * tmask(:,:,:)
178	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
179	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
180	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
181	!
182	END SUBROUTINE istate_t_s
183
184
185	SUBROUTINE istate_eel
186	!!----------------------------------------------------------------------
187	!! * ROUTINE istate_eel *
188	!!
189	!! ** Purpose : Initialization of the dynamics and tracers for EEL R5
190	!! configuration (channel with or without a topographic bump)
191	!!
192	!! ** Method : - set temprature field
193	!! - set salinity field
194	!! - set velocity field including horizontal divergence
195	!! and relative vorticity fields
196	!!----------------------------------------------------------------------
197	USE divcur ! hor. divergence & rel. vorticity (div_cur routine)
198	USE iom
199
200	INTEGER :: inum ! temporary logical unit
201	INTEGER :: ji, jj, jk ! dummy loop indices
202	INTEGER :: ijloc
203	REAL(wp) :: zh1, zh2, zslope, zcst, zfcor ! temporary scalars
204	REAL(wp) :: zt1 = 15._wp ! surface temperature value (EEL R5)
205	REAL(wp) :: zt2 = 5._wp ! bottom temperature value (EEL R5)
206	REAL(wp) :: zsal = 35.0_wp ! constant salinity (EEL R2, R5 and R6)
207	REAL(wp) :: zueel = 0.1_wp ! constant uniform zonal velocity (EEL R5)
208	REAL(wp), DIMENSION(jpiglo,jpjglo) :: zssh ! initial ssh over the global domain
209	!!----------------------------------------------------------------------
210
211	SELECT CASE ( jp_cfg )
212	! ! ====================
213	CASE ( 5 ) ! EEL R5 configuration
214	! ! ====================
215	!
216	! set temperature field with a linear profile
217	! -------------------------------------------
218	IF(lwp) WRITE(numout,*)
219	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile'
220	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
221	!
222	zh1 = gdept_0( 1 )
223	zh2 = gdept_0(jpkm1)
224	!
225	zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )
226	zcst = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 )
227	!
228	DO jk = 1, jpk
229	tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
230	tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
231	END DO
232	!
233	IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &
234	& 1 , jpi , 5 , 1 , jpk , &
235	& 1 , 1. , numout )
236	!
237	! set salinity field to a constant value
238	! --------------------------------------
239	IF(lwp) WRITE(numout,*)
240	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
241	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
242	!
243	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
244	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
245	!
246	! set the dynamics: U,V, hdiv, rot (and ssh if necessary)
247	! ----------------
248	! Start EEL5 configuration with barotropic geostrophic velocities
249	! according the sshb and sshn SSH imposed.
250	! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)
251	! we use the Coriolis frequency at mid-channel.
252	ub(:,:,:) = zueel * umask(:,:,:)
253	un(:,:,:) = ub(:,:,:)
254	ijloc = mj0(INT(jpjglo-1)/2)
255	zfcor = ff(1,ijloc)
256	!
257	DO jj = 1, jpjglo
258	zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)zueele1t(1,1)*zfcor/grav
259	END DO
260	!
261	IF(lwp) THEN
262	WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel
263	WRITE(numout,*) ' Geostrophic SSH profile as a function of y:'
264	WRITE(numout,'(12(1x,f6.2))') zssh(1,:)
265	ENDIF
266	!
267	DO jj = 1, nlcj
268	DO ji = 1, nlci
269	sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1)
270	END DO
271	END DO
272	sshb(nlci+1:jpi, : ) = 0.e0 ! set to zero extra mpp columns
273	sshb( : ,nlcj+1:jpj) = 0.e0 ! set to zero extra mpp rows
274	!
275	sshn(:,:) = sshb(:,:) ! set now ssh to the before value
276	!
277	IF( nn_rstssh /= 0 ) THEN
278	nn_rstssh = 0 ! hand-made initilization of ssh
279	CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )
280	ENDIF
281	!
282	CALL div_cur( nit000 ) ! horizontal divergence and relative vorticity (curl)
283	! N.B. the vertical velocity will be computed from the horizontal divergence field
284	! in istate by a call to wzv routine
285
286
287	! ! ==========================
288	CASE ( 2 , 6 ) ! EEL R2 or R6 configuration
289	! ! ==========================
290	!
291	! set temperature field with a NetCDF file
292	! ----------------------------------------
293	IF(lwp) WRITE(numout,*)
294	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file'
295	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
296	!
297	CALL iom_open ( 'eel.initemp', inum )
298	CALL iom_get ( inum, jpdom_data, 'initemp', tsb(:,:,:,jp_tem) ) ! read before temprature (tb)
299	CALL iom_close( inum )
300	!
301	tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) ! set nox temperature to tb
302	!
303	IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &
304	& 1 , jpi , 5 , 1 , jpk , &
305	& 1 , 1. , numout )
306	!
307	! set salinity field to a constant value
308	! --------------------------------------
309	IF(lwp) WRITE(numout,*)
310	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
311	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
312	!
313	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
314	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
315	!
316	! ! ===========================
317	CASE DEFAULT ! NONE existing configuration
318	! ! ===========================
319	WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded'
320	CALL ctl_stop( ctmp1 )
321	!
322	END SELECT
323	!
324	END SUBROUTINE istate_eel
325
326
327	SUBROUTINE istate_gyre
328	!!----------------------------------------------------------------------
329	!! * ROUTINE istate_gyre *
330	!!
331	!! ** Purpose : Initialization of the dynamics and tracers for GYRE
332	!! configuration (double gyre with rotated domain)
333	!!
334	!! ** Method : - set temprature field
335	!! - set salinity field
336	!!----------------------------------------------------------------------
337	INTEGER :: ji, jj, jk ! dummy loop indices
338	INTEGER :: inum ! temporary logical unit
339	INTEGER, PARAMETER :: ntsinit = 0 ! (0/1) (analytical/input data files) T&S initialization
340	!!----------------------------------------------------------------------
341
342	SELECT CASE ( ntsinit)
343
344	CASE ( 0 ) ! analytical T/S profil deduced from LEVITUS
345	IF(lwp) WRITE(numout,*)
346	IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS '
347	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
348
349	DO jk = 1, jpk
350	DO jj = 1, jpj
351	DO ji = 1, jpi
352	tsn(ji,jj,jk,jp_tem) = ( 16. - 12. * TANH( (fsdept(ji,jj,jk) - 400) / 700 ) ) &
353	& * (-TANH( (500-fsdept(ji,jj,jk)) / 150 ) + 1) / 2 &
354	& + ( 15. * ( 1. - TANH( (fsdept(ji,jj,jk)-50.) / 1500.) ) &
355	& - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.) &
356	& + 7. * (1500. - fsdept(ji,jj,jk)) / 1500. ) &
357	& * (-TANH( (fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
358	tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
359	tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem)
360
361	tsn(ji,jj,jk,jp_sal) = ( 36.25 - 1.13 * TANH( (fsdept(ji,jj,jk) - 305) / 460 ) ) &
362	& * (-TANH((500 - fsdept(ji,jj,jk)) / 150) + 1) / 2 &
363	& + ( 35.55 + 1.25 * (5000. - fsdept(ji,jj,jk)) / 5000. &
364	& - 1.62 * TANH( (fsdept(ji,jj,jk) - 60. ) / 650. ) &
365	& + 0.2 * TANH( (fsdept(ji,jj,jk) - 35. ) / 100. ) &
366	& + 0.2 * TANH( (fsdept(ji,jj,jk) - 1000.) / 5000.) ) &
367	& * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
368	tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
369	tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)
370	END DO
371	END DO
372	END DO
373
374	CASE ( 1 ) ! T/S data fields read in dta_tem.nc/data_sal.nc files
375	IF(lwp) WRITE(numout,*)
376	IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files'
377	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
378	IF(lwp) WRITE(numout,*) ' NetCDF FORMAT'
379
380	! Read temperature field
381	! ----------------------
382	CALL iom_open ( 'data_tem', inum )
383	CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) )
384	CALL iom_close( inum )
385
386	tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
387	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
388
389	! Read salinity field
390	! -------------------
391	CALL iom_open ( 'data_sal', inum )
392	CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) )
393	CALL iom_close( inum )
394
395	tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
396	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
397
398	END SELECT
399
400	IF(lwp) THEN
401	WRITE(numout,*)
402	WRITE(numout,*) ' Initial temperature and salinity profiles:'
403	WRITE(numout, "(9x,' level gdept_0 temperature salinity ')" )
404	WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_0(jk), tsn(2,2,jk,jp_tem), tsn(2,2,jk,jp_sal), jk = 1, jpk )
405	ENDIF
406
407	END SUBROUTINE istate_gyre
408
409
410	SUBROUTINE istate_uvg
411	!!----------------------------------------------------------------------
412	!! * ROUTINE istate_uvg *
413	!!
414	!! ** Purpose : Compute the geostrophic velocities from (tn,sn) fields
415	!!
416	!! ** Method : Using the hydrostatic hypothesis the now hydrostatic
417	!! pressure is computed by integrating the in-situ density from the
418	!! surface to the bottom.
419	!! p=integral [ rau*g dz ]
420	!!----------------------------------------------------------------------
421	USE dynspg ! surface pressure gradient (dyn_spg routine)
422	USE divcur ! hor. divergence & rel. vorticity (div_cur routine)
423	USE lbclnk ! ocean lateral boundary condition (or mpp link)
424
425	INTEGER :: ji, jj, jk ! dummy loop indices
426	INTEGER :: indic ! ???
427	REAL(wp) :: zmsv, zphv, zmsu, zphu, zalfg ! temporary scalars
428	REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn
429	!!----------------------------------------------------------------------
430	!
431	CALL wrk_alloc( jpi, jpj, jpk, zprn)
432	!
433	IF(lwp) WRITE(numout,*)
434	IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'
435	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
436
437	! Compute the now hydrostatic pressure
438	! ------------------------------------
439
440	zalfg = 0.5 * grav * rau0
441
442	zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value
443
444	DO jk = 2, jpkm1 ! Vertical integration from the surface
445	zprn(:,:,jk) = zprn(:,:,jk-1) &
446	& + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
447	END DO
448
449	! Compute geostrophic balance
450	! ---------------------------
451	DO jk = 1, jpkm1
452	DO jj = 2, jpjm1
453	DO ji = fs_2, fs_jpim1 ! vertor opt.
454	zmsv = 1. / MAX( umask(ji-1,jj+1,jk) + umask(ji ,jj+1,jk) &
455	+ umask(ji-1,jj ,jk) + umask(ji ,jj ,jk) , 1. )
456	zphv = ( zprn(ji ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1) &
457	+ ( zprn(ji+1,jj+1,jk) - zprn(ji ,jj+1,jk) ) * umask(ji ,jj+1,jk) / e1u(ji ,jj+1) &
458	+ ( zprn(ji ,jj ,jk) - zprn(ji-1,jj ,jk) ) * umask(ji-1,jj ,jk) / e1u(ji-1,jj ) &
459	+ ( zprn(ji+1,jj ,jk) - zprn(ji ,jj ,jk) ) * umask(ji ,jj ,jk) / e1u(ji ,jj )
460	zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk)
461
462	zmsu = 1. / MAX( vmask(ji+1,jj ,jk) + vmask(ji ,jj ,jk) &
463	+ vmask(ji+1,jj-1,jk) + vmask(ji ,jj-1,jk) , 1. )
464	zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj ,jk) ) * vmask(ji+1,jj ,jk) / e2v(ji+1,jj ) &
465	+ ( zprn(ji ,jj+1,jk) - zprn(ji ,jj ,jk) ) * vmask(ji ,jj ,jk) / e2v(ji ,jj ) &
466	+ ( zprn(ji+1,jj ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1) &
467	+ ( zprn(ji ,jj ,jk) - zprn(ji ,jj-1,jk) ) * vmask(ji ,jj-1,jk) / e2v(ji ,jj-1)
468	zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk)
469
470	! Compute the geostrophic velocities
471	un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji ,jj-1) )
472	vn(ji,jj,jk) = 2. * zphv / ( ff(ji,jj) + ff(ji-1,jj ) )
473	END DO
474	END DO
475	END DO
476
477	IF(lwp) WRITE(numout,*) ' we force to zero bottom velocity'
478
479	! Susbtract the bottom velocity (level jpk-1 for flat bottom case)
480	! to have a zero bottom velocity
481
482	DO jk = 1, jpkm1
483	un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk)
484	vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk)
485	END DO
486
487	CALL lbc_lnk( un, 'U', -1. )
488	CALL lbc_lnk( vn, 'V', -1. )
489
490	ub(:,:,:) = un(:,:,:)
491	vb(:,:,:) = vn(:,:,:)
492
493	! WARNING !!!!!
494	! after initializing u and v, we need to calculate the initial streamfunction bsf.
495	! Otherwise, only the trend will be computed and the model will blow up (inconsistency).
496	! to do that, we call dyn_spg with a special trick:
497	! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the
498	! right value assuming the velocities have been set up in one time step.
499	! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)
500	! sets up s false trend to calculate the barotropic streamfunction.
501
502	ua(:,:,:) = ub(:,:,:) / rdt
503	va(:,:,:) = vb(:,:,:) / rdt
504
505	! calls dyn_spg. we assume euler time step, starting from rest.
506	indic = 0
507	CALL dyn_spg( nit000, indic ) ! surface pressure gradient
508
509	! the new velocity is ua*rdt
510
511	CALL lbc_lnk( ua, 'U', -1. )
512	CALL lbc_lnk( va, 'V', -1. )
513
514	ub(:,:,:) = ua(:,:,:) * rdt
515	vb(:,:,:) = va(:,:,:) * rdt
516	ua(:,:,:) = 0.e0
517	va(:,:,:) = 0.e0
518	un(:,:,:) = ub(:,:,:)
519	vn(:,:,:) = vb(:,:,:)
520
521	! Compute the divergence and curl
522
523	CALL div_cur( nit000 ) ! now horizontal divergence and curl
524
525	hdivb(:,:,:) = hdivn(:,:,:) ! set the before to the now value
526	rotb (:,:,:) = rotn (:,:,:) ! set the before to the now value
527	!
528	CALL wrk_dealloc( jpi, jpj, jpk, zprn)
529	!
530	END SUBROUTINE istate_uvg
531
532	!!=====================================================================
533	END MODULE istate

Download in other formats:

Original Format

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

Context Navigation

Ticket #1126: istate.F90

Download in other formats: