Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

istate.F90 @ 8011

Last change on this file since 8011 was 8011, checked in by jchanut, 7 years ago
VORTEX test case for AGRIF
File size: 31.3 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	!! 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 c1d ! 1D vertical configuration
29	USE daymod ! calendar
30	USE eosbn2 ! eq. of state, Brunt Vaisala frequency (eos routine)
31	USE ldftra_oce ! ocean active tracers: lateral physics
32	USE zdf_oce ! ocean vertical physics
33	USE phycst ! physical constants
34	USE dtatsd ! data temperature and salinity (dta_tsd routine)
35	USE dtauvd ! data: U & V current (dta_uvd routine)
36	USE zpshde ! partial step: hor. derivative (zps_hde routine)
37	USE eosbn2 ! equation of state (eos bn2 routine)
38	USE domvvl ! varying vertical mesh
39	USE dynspg_oce ! pressure gradient schemes
40	USE dynspg_flt ! filtered free surface
41	USE sol_oce ! ocean solver variables
42	!
43	USE in_out_manager ! I/O manager
44	USE iom ! I/O library
45	USE lib_mpp ! MPP library
46	USE restart ! restart
47	USE wrk_nemo ! Memory allocation
48	USE timing ! Timing
49
50	IMPLICIT NONE
51	PRIVATE
52
53	PUBLIC istate_init ! routine called by step.F90
54
55	!! * Substitutions
56	# include "domzgr_substitute.h90"
57	# include "vectopt_loop_substitute.h90"
58	!!----------------------------------------------------------------------
59	!! NEMO/OPA 3.7 , NEMO Consortium (2014)
60	!! $Id: istate.F90 5332 2015-06-01 16:59:09Z mathiot $
61	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
62	!!----------------------------------------------------------------------
63	CONTAINS
64
65	SUBROUTINE istate_init
66	!!----------------------------------------------------------------------
67	!! * ROUTINE istate_init *
68	!!
69	!! ** Purpose : Initialization of the dynamics and tracer fields.
70	!!----------------------------------------------------------------------
71	INTEGER :: ji, jj, jk ! dummy loop indices
72	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zuvd ! U & V data workspace
73	!!----------------------------------------------------------------------
74	!
75	IF( nn_timing == 1 ) CALL timing_start('istate_init')
76	!
77
78	IF(lwp) WRITE(numout,*) ' '
79	IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
80	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
81
82	CALL dta_tsd_init ! Initialisation of T & S input data
83	IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data
84
85	rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk
86	rn2b (:,:,: ) = 0._wp ; rn2 (:,:,: ) = 0._wp ! set one for all to 0 at levels 1 and jpk
87	tsa (:,:,:,:) = 0._wp ! set one for all to 0 at level jpk
88	rab_b(:,:,:,:) = 0._wp ; rab_n(:,:,:,:) = 0._wp ! set one for all to 0 at level jpk
89
90	IF( ln_rstart ) THEN ! Restart from a file
91	! ! -------------------
92	CALL rst_read ! Read the restart file
93	CALL day_init ! model calendar (using both namelist and restart infos)
94	ELSE
95	! ! Start from rest
96	! ! ---------------
97	numror = 0 ! define numror = 0 -> no restart file to read
98	neuler = 0 ! Set time-step indicator at nit000 (euler forward)
99	CALL day_init ! model calendar (using both namelist and restart infos)
100	! ! Initialization of ocean to zero
101	! before fields ! now fields
102	sshb (:,:) = 0._wp ; sshn (:,:) = 0._wp
103	ub (:,:,:) = 0._wp ; un (:,:,:) = 0._wp
104	vb (:,:,:) = 0._wp ; vn (:,:,:) = 0._wp
105	rotb (:,:,:) = 0._wp ; rotn (:,:,:) = 0._wp
106	hdivb(:,:,:) = 0._wp ; hdivn(:,:,:) = 0._wp
107	!
108	IF( cp_cfg == 'eel' ) THEN
109	CALL istate_eel ! EEL configuration : start from pre-defined U,V T-S fields
110	ELSEIF( cp_cfg == 'gyre' ) THEN
111	CALL istate_gyre ! GYRE configuration : start from pre-defined T-S fields
112	ELSEIF( cp_cfg == 'vortex' ) THEN
113	CALL istate_vortex ! GYRE configuration : start from pre-defined T-S fields
114	ELSE ! Initial T-S, U-V fields read in files
115	IF ( ln_tsd_init ) THEN ! read 3D T and S data at nit000
116	CALL dta_tsd( nit000, tsb )
117	tsn(:,:,:,:) = tsb(:,:,:,:)
118	!
119	ELSE ! Initial T-S fields defined analytically
120	CALL istate_t_s
121	ENDIF
122	IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
123	CALL wrk_alloc( jpi, jpj, jpk, 2, zuvd )
124	CALL dta_uvd( nit000, zuvd )
125	ub(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:)
126	vb(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:)
127	CALL wrk_dealloc( jpi, jpj, jpk, 2, zuvd )
128	ENDIF
129	ENDIF
130	!
131	CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! before potential and in situ densities
132	#if ! defined key_c1d
133	IF( ln_zps .AND. .NOT. ln_isfcav) &
134	& CALL zps_hde ( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient
135	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
136	IF( ln_zps .AND. ln_isfcav) &
137	& CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF)
138	& rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , &
139	& gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level
140	#endif
141	!
142	! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
143	IF( lk_vvl ) THEN
144	DO jk = 1, jpk
145	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
146	ENDDO
147	ENDIF
148	!
149	ENDIF
150	!
151	IF( lk_agrif ) THEN ! read free surface arrays in restart file
152	IF( ln_rstart ) THEN
153	IF( lk_dynspg_flt ) THEN ! read or initialize the following fields
154	! ! gcx, gcxb for agrif_opa_init
155	IF( sol_oce_alloc() > 0 ) CALL ctl_stop('agrif sol_oce_alloc: allocation of arrays failed')
156	CALL flt_rst( nit000, 'READ' )
157	ENDIF
158	ENDIF ! explicit case not coded yet with AGRIF
159	ENDIF
160	!
161	!
162	! Initialize "now" and "before" barotropic velocities:
163	! Do it whatever the free surface method, these arrays
164	! being eventually used
165	!
166	!
167	un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
168	ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
169	!
170	DO jk = 1, jpkm1
171	DO jj = 1, jpj
172	DO ji = 1, jpi
173	un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
174	vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
175	!
176	ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
177	vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
178	END DO
179	END DO
180	END DO
181	!
182	un_b(:,:) = un_b(:,:) * hur (:,:)
183	vn_b(:,:) = vn_b(:,:) * hvr (:,:)
184	!
185	ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
186	vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
187	!
188	!
189	IF( nn_timing == 1 ) CALL timing_stop('istate_init')
190	!
191	END SUBROUTINE istate_init
192
193
194	SUBROUTINE istate_t_s
195	!!---------------------------------------------------------------------
196	!! * ROUTINE istate_t_s *
197	!!
198	!! ** Purpose : Intialization of the temperature field with an
199	!! analytical profile or a file (i.e. in EEL configuration)
200	!!
201	!! ** Method : - temperature: use Philander analytic profile
202	!! - salinity : use to a constant value 35.5
203	!!
204	!! References : Philander ???
205	!!----------------------------------------------------------------------
206	INTEGER :: ji, jj, jk
207	REAL(wp) :: zsal = 35.50
208	!!----------------------------------------------------------------------
209	!
210	IF(lwp) WRITE(numout,*)
211	IF(lwp) WRITE(numout,*) 'istate_t_s : Philander s initial temperature profile'
212	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ and constant salinity (',zsal,' psu)'
213	!
214	DO jk = 1, jpk
215	tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) ) &
216	& + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.) ) * tmask(:,:,jk)
217	tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
218	END DO
219	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
220	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
221	!
222	END SUBROUTINE istate_t_s
223
224
225	SUBROUTINE istate_eel
226	!!----------------------------------------------------------------------
227	!! * ROUTINE istate_eel *
228	!!
229	!! ** Purpose : Initialization of the dynamics and tracers for EEL R5
230	!! configuration (channel with or without a topographic bump)
231	!!
232	!! ** Method : - set temprature field
233	!! - set salinity field
234	!! - set velocity field including horizontal divergence
235	!! and relative vorticity fields
236	!!----------------------------------------------------------------------
237	USE divcur ! hor. divergence & rel. vorticity (div_cur routine)
238	USE iom
239	!
240	INTEGER :: inum ! temporary logical unit
241	INTEGER :: ji, jj, jk ! dummy loop indices
242	INTEGER :: ijloc
243	REAL(wp) :: zh1, zh2, zslope, zcst, zfcor ! temporary scalars
244	REAL(wp) :: zt1 = 15._wp ! surface temperature value (EEL R5)
245	REAL(wp) :: zt2 = 5._wp ! bottom temperature value (EEL R5)
246	REAL(wp) :: zsal = 35.0_wp ! constant salinity (EEL R2, R5 and R6)
247	REAL(wp) :: zueel = 0.1_wp ! constant uniform zonal velocity (EEL R5)
248	REAL(wp), DIMENSION(jpiglo,jpjglo) :: zssh ! initial ssh over the global domain
249	!!----------------------------------------------------------------------
250	!
251	SELECT CASE ( jp_cfg )
252	! ! ====================
253	CASE ( 5 ) ! EEL R5 configuration
254	! ! ====================
255	!
256	! set temperature field with a linear profile
257	! -------------------------------------------
258	IF(lwp) WRITE(numout,*)
259	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile'
260	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
261	!
262	zh1 = gdept_1d( 1 )
263	zh2 = gdept_1d(jpkm1)
264	!
265	zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )
266	zcst = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 )
267	!
268	DO jk = 1, jpk
269	tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
270	tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
271	END DO
272	!
273	IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &
274	& 1 , jpi , 5 , 1 , jpk , &
275	& 1 , 1. , numout )
276	!
277	! set salinity field to a constant value
278	! --------------------------------------
279	IF(lwp) WRITE(numout,*)
280	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
281	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
282	!
283	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
284	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
285	!
286	! set the dynamics: U,V, hdiv, rot (and ssh if necessary)
287	! ----------------
288	! Start EEL5 configuration with barotropic geostrophic velocities
289	! according the sshb and sshn SSH imposed.
290	! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)
291	! we use the Coriolis frequency at mid-channel.
292	ub(:,:,:) = zueel * umask(:,:,:)
293	un(:,:,:) = ub(:,:,:)
294	ijloc = mj0(INT(jpjglo-1)/2)
295	zfcor = ff(1,ijloc)
296	!
297	DO jj = 1, jpjglo
298	zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)zueele1t(1,1)*zfcor/grav
299	END DO
300	!
301	IF(lwp) THEN
302	WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel
303	WRITE(numout,*) ' Geostrophic SSH profile as a function of y:'
304	WRITE(numout,'(12(1x,f6.2))') zssh(1,:)
305	ENDIF
306	!
307	DO jj = 1, nlcj
308	DO ji = 1, nlci
309	sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1)
310	END DO
311	END DO
312	sshb(nlci+1:jpi, : ) = 0.e0 ! set to zero extra mpp columns
313	sshb( : ,nlcj+1:jpj) = 0.e0 ! set to zero extra mpp rows
314	!
315	sshn(:,:) = sshb(:,:) ! set now ssh to the before value
316	!
317	IF( nn_rstssh /= 0 ) THEN
318	nn_rstssh = 0 ! hand-made initilization of ssh
319	CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )
320	ENDIF
321	!
322	CALL div_cur( nit000 ) ! horizontal divergence and relative vorticity (curl)
323	! N.B. the vertical velocity will be computed from the horizontal divergence field
324	! in istate by a call to wzv routine
325
326
327	! ! ==========================
328	CASE ( 2 , 6 ) ! EEL R2 or R6 configuration
329	! ! ==========================
330	!
331	! set temperature field with a NetCDF file
332	! ----------------------------------------
333	IF(lwp) WRITE(numout,*)
334	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file'
335	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
336	!
337	CALL iom_open ( 'eel.initemp', inum )
338	CALL iom_get ( inum, jpdom_data, 'initemp', tsb(:,:,:,jp_tem) ) ! read before temprature (tb)
339	CALL iom_close( inum )
340	!
341	tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) ! set nox temperature to tb
342	!
343	IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi , jpj , jpk , jpj/2 , &
344	& 1 , jpi , 5 , 1 , jpk , &
345	& 1 , 1. , numout )
346	!
347	! set salinity field to a constant value
348	! --------------------------------------
349	IF(lwp) WRITE(numout,*)
350	IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
351	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
352	!
353	tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
354	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
355	!
356	! ! ===========================
357	CASE DEFAULT ! NONE existing configuration
358	! ! ===========================
359	WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded'
360	CALL ctl_stop( ctmp1 )
361	!
362	END SELECT
363	!
364	END SUBROUTINE istate_eel
365
366
367	SUBROUTINE istate_gyre
368	!!----------------------------------------------------------------------
369	!! * ROUTINE istate_gyre *
370	!!
371	!! ** Purpose : Initialization of the dynamics and tracers for GYRE
372	!! configuration (double gyre with rotated domain)
373	!!
374	!! ** Method : - set temprature field
375	!! - set salinity field
376	!!----------------------------------------------------------------------
377	INTEGER :: ji, jj, jk ! dummy loop indices
378	INTEGER :: inum ! temporary logical unit
379	INTEGER, PARAMETER :: ntsinit = 0 ! (0/1) (analytical/input data files) T&S initialization
380	!!----------------------------------------------------------------------
381	!
382	SELECT CASE ( ntsinit)
383	!
384	CASE ( 0 ) ! analytical T/S profil deduced from LEVITUS
385	IF(lwp) WRITE(numout,*)
386	IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS '
387	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
388	!
389	DO jk = 1, jpk
390	DO jj = 1, jpj
391	DO ji = 1, jpi
392	tsn(ji,jj,jk,jp_tem) = ( 16. - 12. * TANH( (fsdept(ji,jj,jk) - 400) / 700 ) ) &
393	& * (-TANH( (500-fsdept(ji,jj,jk)) / 150 ) + 1) / 2 &
394	& + ( 15. * ( 1. - TANH( (fsdept(ji,jj,jk)-50.) / 1500.) ) &
395	& - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.) &
396	& + 7. * (1500. - fsdept(ji,jj,jk)) / 1500. ) &
397	& * (-TANH( (fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
398	tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
399	tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem)
400
401	tsn(ji,jj,jk,jp_sal) = ( 36.25 - 1.13 * TANH( (fsdept(ji,jj,jk) - 305) / 460 ) ) &
402	& * (-TANH((500 - fsdept(ji,jj,jk)) / 150) + 1) / 2 &
403	& + ( 35.55 + 1.25 * (5000. - fsdept(ji,jj,jk)) / 5000. &
404	& - 1.62 * TANH( (fsdept(ji,jj,jk) - 60. ) / 650. ) &
405	& + 0.2 * TANH( (fsdept(ji,jj,jk) - 35. ) / 100. ) &
406	& + 0.2 * TANH( (fsdept(ji,jj,jk) - 1000.) / 5000.) ) &
407	& * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
408	tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
409	tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)
410	END DO
411	END DO
412	END DO
413	!
414	CASE ( 1 ) ! T/S data fields read in dta_tem.nc/data_sal.nc files
415	IF(lwp) WRITE(numout,*)
416	IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files'
417	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
418	IF(lwp) WRITE(numout,*) ' NetCDF FORMAT'
419
420	! Read temperature field
421	! ----------------------
422	CALL iom_open ( 'data_tem', inum )
423	CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) )
424	CALL iom_close( inum )
425
426	tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
427	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
428
429	! Read salinity field
430	! -------------------
431	CALL iom_open ( 'data_sal', inum )
432	CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) )
433	CALL iom_close( inum )
434
435	tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
436	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
437	!
438	END SELECT
439	!
440	IF(lwp) THEN
441	WRITE(numout,*)
442	WRITE(numout,*) ' Initial temperature and salinity profiles:'
443	WRITE(numout, "(9x,' level gdept_1d temperature salinity ')" )
444	WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_1d(jk), tsn(2,2,jk,jp_tem), tsn(2,2,jk,jp_sal), jk = 1, jpk )
445	ENDIF
446	!
447	END SUBROUTINE istate_gyre
448
449	SUBROUTINE istate_vortex
450	!!----------------------------------------------------------------------
451	!! * ROUTINE istate_vortex *
452	!!
453	!! ** Purpose : Initialization of the dynamics and tracers for
454	!! baroclinic vortex
455	!!
456	!! ** Method : - set ssh
457	!! - set temperature/salinity fields
458	!! - set velocities
459	!!
460	!! ** Ref : Penven et al. (2006) : Evaluation and application of
461	!! the ROMS 1-way embedding procedure to the central
462	!! california upwelling system. Ocean Mod., 12, 1–2, 157–187.
463	!!----------------------------------------------------------------------
464	INTEGER :: ji, jj, jk ! dummy loop indices
465	REAL(wp) :: zx, zy, zP0, zumax, zlambda, zn2, zf0, zH, zrho1, za, zf
466	REAL(wp) :: zdt, zdu, zdv
467	!!----------------------------------------------------------------------
468	!
469	!
470	zf0 = 2._wp * omega * SIN( rad * ppgphi0 )
471	zumax = 1._wp ! Anticyclonic: set zumax=-1 for cyclonic
472	zlambda = SQRT(2._wp)*60.e3
473	zn2 = 3.e-3**2
474	zH = 0.5_wp * rn_bathy
475	!
476	zP0 = rau0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)
477	! Sea level:
478	za = -zP0 * (1._wp-EXP(-zH)) / (grav*(zH-1._wp + EXP(-zH)))
479	DO ji=1, jpi
480	DO jj=1, jpj
481	zx = glamt(ji,jj) * 1.e3
482	zy = gphit(ji,jj) * 1.e3
483	zrho1 = rau0 + za * EXP(-(zx2+zy2)/zlambda**2)
484	sshn(ji,jj) = zP0 * EXP(-(zx2+zy2)/zlambda*2)/(zrho1grav) * tmask(ji,jj,1)
485	END DO
486	END DO
487	sshb(:,:) = sshn(:,:)
488	!
489	! temperature:
490	DO ji=1, jpi
491	DO jj=1, jpj
492	zx = glamt(ji,jj) * 1.e3
493	zy = gphit(ji,jj) * 1.e3
494	DO jk=1,jpk
495	zdt = fsdept(ji,jj,jk)
496	zrho1 = rau0 * (1._wp + zn2*zdt/grav)
497	IF (zdt < zH) THEN
498	zrho1 = zrho1 - zP0 * (1._wp-EXP(zdt-zH)) &
499	& * EXP(-(zx2+zy2)/zlambda*2) / (grav(zH -1._wp + exp(-zH)));
500	ENDIF
501	tsn(ji,jj,jk,jp_tem) = (20._wp + (rau0-zrho1) / 0.28_wp) * tmask(ji,jj,jk)
502	END DO
503	END DO
504	END DO
505	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
506	!
507	! salinity:
508	tsn(:,:,:,jp_sal) = 35._wp
509	tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
510	!
511	! velocities:
512	za = 2._wp * zP0 / (zf0 * rau0 * zlambda**2)
513	DO ji=1, jpim1
514	DO jj=1, jpj
515	zx = glamu(ji,jj) * 1.e3
516	zy = gphiu(ji,jj) * 1.e3
517	DO jk=1, jpk
518	zdu = 0.5_wp * (fsdept(ji ,jj,jk) + fsdept(ji+1,jj,jk))
519	zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH))
520	IF (zdu < zH) THEN
521	un(ji,jj,jk) = (za * zf * zy * EXP(-(zx2+zy2)/zlambda*2)) umask(ji,jj,jk)
522	ELSE
523	un(ji,jj,jk) = 0._wp
524	ENDIF
525	END DO
526	END DO
527	END DO
528	!
529	DO ji=1, jpi
530	DO jj=1, jpjm1
531	zx = glamv(ji,jj) * 1.e3
532	zy = gphiv(ji,jj) * 1.e3
533	DO jk=1, jpk
534	zdv = 0.5_wp * (fsdept(ji ,jj,jk) + fsdept(ji,jj+1,jk))
535	zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH))
536	IF (zdv < zH) THEN
537	vn(ji,jj,jk) = -(za * zf * zx * EXP(-(zx2+zy2)/zlambda*2)) umask(ji,jj,jk)
538	ELSE
539	vn(ji,jj,jk) = 0._wp
540	ENDIF
541	END DO
542	END DO
543	END DO
544
545	CALL lbc_lnk( un, 'U', -1. )
546	CALL lbc_lnk( vn, 'V', -1. )
547
548	ub(:,:,:) = un(:,:,:)
549	vb(:,:,:) = vn(:,:,:)
550
551	END SUBROUTINE istate_vortex
552
553	SUBROUTINE istate_uvg
554	!!----------------------------------------------------------------------
555	!! * ROUTINE istate_uvg *
556	!!
557	!! ** Purpose : Compute the geostrophic velocities from (tn,sn) fields
558	!!
559	!! ** Method : Using the hydrostatic hypothesis the now hydrostatic
560	!! pressure is computed by integrating the in-situ density from the
561	!! surface to the bottom.
562	!! p=integral [ rau*g dz ]
563	!!----------------------------------------------------------------------
564	USE dynspg ! surface pressure gradient (dyn_spg routine)
565	USE divcur ! hor. divergence & rel. vorticity (div_cur routine)
566	USE lbclnk ! ocean lateral boundary condition (or mpp link)
567	!
568	INTEGER :: ji, jj, jk ! dummy loop indices
569	INTEGER :: indic ! ???
570	REAL(wp) :: zmsv, zphv, zmsu, zphu, zalfg ! temporary scalars
571	REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn
572	!!----------------------------------------------------------------------
573	!
574	CALL wrk_alloc( jpi, jpj, jpk, zprn)
575	!
576	IF(lwp) WRITE(numout,*)
577	IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'
578	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
579
580	! Compute the now hydrostatic pressure
581	! ------------------------------------
582
583	zalfg = 0.5 * grav * rau0
584
585	zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value
586
587	DO jk = 2, jpkm1 ! Vertical integration from the surface
588	zprn(:,:,jk) = zprn(:,:,jk-1) &
589	& + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
590	END DO
591
592	! Compute geostrophic balance
593	! ---------------------------
594	DO jk = 1, jpkm1
595	DO jj = 2, jpjm1
596	DO ji = fs_2, fs_jpim1 ! vertor opt.
597	zmsv = 1. / MAX( umask(ji-1,jj+1,jk) + umask(ji ,jj+1,jk) &
598	+ umask(ji-1,jj ,jk) + umask(ji ,jj ,jk) , 1. )
599	zphv = ( zprn(ji ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1) &
600	+ ( zprn(ji+1,jj+1,jk) - zprn(ji ,jj+1,jk) ) * umask(ji ,jj+1,jk) / e1u(ji ,jj+1) &
601	+ ( zprn(ji ,jj ,jk) - zprn(ji-1,jj ,jk) ) * umask(ji-1,jj ,jk) / e1u(ji-1,jj ) &
602	+ ( zprn(ji+1,jj ,jk) - zprn(ji ,jj ,jk) ) * umask(ji ,jj ,jk) / e1u(ji ,jj )
603	zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk)
604
605	zmsu = 1. / MAX( vmask(ji+1,jj ,jk) + vmask(ji ,jj ,jk) &
606	+ vmask(ji+1,jj-1,jk) + vmask(ji ,jj-1,jk) , 1. )
607	zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj ,jk) ) * vmask(ji+1,jj ,jk) / e2v(ji+1,jj ) &
608	+ ( zprn(ji ,jj+1,jk) - zprn(ji ,jj ,jk) ) * vmask(ji ,jj ,jk) / e2v(ji ,jj ) &
609	+ ( zprn(ji+1,jj ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1) &
610	+ ( zprn(ji ,jj ,jk) - zprn(ji ,jj-1,jk) ) * vmask(ji ,jj-1,jk) / e2v(ji ,jj-1)
611	zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk)
612
613	! Compute the geostrophic velocities
614	un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji ,jj-1) )
615	vn(ji,jj,jk) = 2. * zphv / ( ff(ji,jj) + ff(ji-1,jj ) )
616	END DO
617	END DO
618	END DO
619
620	IF(lwp) WRITE(numout,*) ' we force to zero bottom velocity'
621
622	! Susbtract the bottom velocity (level jpk-1 for flat bottom case)
623	! to have a zero bottom velocity
624
625	DO jk = 1, jpkm1
626	un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk)
627	vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk)
628	END DO
629
630	CALL lbc_lnk( un, 'U', -1. )
631	CALL lbc_lnk( vn, 'V', -1. )
632
633	ub(:,:,:) = un(:,:,:)
634	vb(:,:,:) = vn(:,:,:)
635
636	! WARNING !!!!!
637	! after initializing u and v, we need to calculate the initial streamfunction bsf.
638	! Otherwise, only the trend will be computed and the model will blow up (inconsistency).
639	! to do that, we call dyn_spg with a special trick:
640	! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the
641	! right value assuming the velocities have been set up in one time step.
642	! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)
643	! sets up s false trend to calculate the barotropic streamfunction.
644
645	ua(:,:,:) = ub(:,:,:) / rdt
646	va(:,:,:) = vb(:,:,:) / rdt
647
648	! calls dyn_spg. we assume euler time step, starting from rest.
649	indic = 0
650	CALL dyn_spg( nit000, indic ) ! surface pressure gradient
651
652	! the new velocity is ua*rdt
653
654	CALL lbc_lnk( ua, 'U', -1. )
655	CALL lbc_lnk( va, 'V', -1. )
656
657	ub(:,:,:) = ua(:,:,:) * rdt
658	vb(:,:,:) = va(:,:,:) * rdt
659	ua(:,:,:) = 0.e0
660	va(:,:,:) = 0.e0
661	un(:,:,:) = ub(:,:,:)
662	vn(:,:,:) = vb(:,:,:)
663
664	! Compute the divergence and curl
665
666	CALL div_cur( nit000 ) ! now horizontal divergence and curl
667
668	hdivb(:,:,:) = hdivn(:,:,:) ! set the before to the now value
669	rotb (:,:,:) = rotn (:,:,:) ! set the before to the now value
670	!
671	CALL wrk_dealloc( jpi, jpj, jpk, zprn)
672	!
673	END SUBROUTINE istate_uvg
674
675	!!=====================================================================
676	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/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/CONFIG/VORTEX/MY_SRC/istate.F90 @ 8011

Download in other formats: