Context Navigation

istate.F90 @ 3304

Last change on this file since 3304 was 2104, checked in by cetlod, 14 years ago
update DEV_r2006_merge_TRA_TRC according to review
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 25.5 KB

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

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

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

Context Navigation

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/istate.F90 @ 3304

Download in other formats: