Context Navigation

Back to Ticket #1438

Ticket #1438: istate.F90

File istate.F90, 26.3 KB (added by Robin_Waldman, 9 years ago)

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

Download in other formats:

Original Format

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

Context Navigation

Ticket #1438: istate.F90

Download in other formats: