Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

flodom.F90 @ 10986

Last change on this file since 10986 was 10986, checked in by andmirek, 5 years ago
GMED 462 add flush
File size: 20.4 KB

Line
1	MODULE flodom
2	!!======================================================================
3	!! * MODULE flodom *
4	!! Ocean floats : domain
5	!!======================================================================
6	!! History : OPA ! 1998-07 (Y.Drillet, CLIPPER) Original code
7	!! NEMO 3.3 ! 2011-09 (C.Bricaud,S.Law-Chune Mercator-Ocean): add ARIANE convention + comsecitc changes
8	!!----------------------------------------------------------------------
9	#if defined key_floats
10	!!----------------------------------------------------------------------
11	!! 'key_floats' float trajectories
12	!!----------------------------------------------------------------------
13	!! flo_dom : initialization of floats
14	!! add_new_floats : add new floats (long/lat/depth)
15	!! add_new_ariane_floats : add new floats with araine convention (i/j/k)
16	!! findmesh : compute index of position
17	!! dstnce : compute distance between face mesh and floats
18	!!----------------------------------------------------------------------
19	USE oce ! ocean dynamics and tracers
20	USE dom_oce ! ocean space and time domain
21	USE flo_oce ! ocean drifting floats
22	USE in_out_manager ! I/O manager
23	USE lib_mpp ! distribued memory computing library
24
25	IMPLICIT NONE
26	PRIVATE
27
28	PUBLIC flo_dom ! routine called by floats.F90
29	PUBLIC flo_dom_alloc ! Routine called in floats.F90
30
31	CHARACTER (len=21) :: clname1 = 'init_float' ! floats initialisation filename
32	CHARACTER (len=21) :: clname2 = 'init_float_ariane' ! ariane floats initialisation filename
33
34
35	INTEGER , ALLOCATABLE, DIMENSION(:) :: iimfl, ijmfl, ikmfl ! index mesh of floats
36	INTEGER , ALLOCATABLE, DIMENSION(:) :: idomfl, ivtest, ihtest ! -
37	REAL(wp), ALLOCATABLE, DIMENSION(:) :: zgifl, zgjfl, zgkfl ! distances in indexes
38
39	!!----------------------------------------------------------------------
40	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
41	!! $Id$
42	!! Software governed by the CeCILL license (see ./LICENSE)
43	!!----------------------------------------------------------------------
44	CONTAINS
45
46	SUBROUTINE flo_dom
47	!! ---------------------------------------------------------------------
48	!! * ROUTINE flo_dom *
49	!!
50	!! ** Purpose : Initialisation of floats
51	!!
52	!! ** Method : We put the floats in the domain with the latitude,
53	!! the longitude (degree) and the depth (m).
54	!!----------------------------------------------------------------------
55	INTEGER :: jfl ! dummy loop
56	INTEGER :: inum ! logical unit for file read
57	!!---------------------------------------------------------------------
58
59	! Initialisation with the geographical position or restart
60
61	IF(lwp) THEN
62	WRITE(numout,*) 'flo_dom : compute initial position of floats'
63	WRITE(numout,*) '~~~~~~~~'
64	WRITE(numout,*) ' jpnfl = ',jpnfl
65	IF(lflush) CALL FLUSH(numout)
66	ENDIF
67
68	!-------------------------!
69	! FLOAT RESTART FILE READ !
70	!-------------------------!
71	IF( ln_rstflo )THEN
72
73	IF(lwp) WRITE(numout,*) ' float restart file read'
74
75	! open the restart file
76	!----------------------
77	CALL ctl_opn( inum, 'restart_float', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
78
79	! read of the restart file
80	READ(inum,*) ( tpifl (jfl), jfl=1, jpnrstflo), &
81	( tpjfl (jfl), jfl=1, jpnrstflo), &
82	( tpkfl (jfl), jfl=1, jpnrstflo), &
83	( nisobfl(jfl), jfl=1, jpnrstflo), &
84	( ngrpfl (jfl), jfl=1, jpnrstflo)
85	CLOSE(inum)
86
87	! if we want a surface drift ( like PROVOR floats )
88	IF( ln_argo ) nisobfl(1:jpnrstflo) = 0
89
90	! It is possible to add new floats.
91	!---------------------------------
92	IF( jpnfl > jpnrstflo )THEN
93
94	IF(lwp) THEN
95	WRITE(numout,*) ' add new floats'
96	IF(lflush) CALL FLUSH(numout)
97	ENDIF
98
99	IF( ln_ariane )THEN !Add new floats with ariane convention
100	CALL flo_add_new_ariane_floats(jpnrstflo+1,jpnfl)
101	ELSE !Add new floats with long/lat convention
102	CALL flo_add_new_floats(jpnrstflo+1,jpnfl)
103	ENDIF
104	ENDIF
105
106	!--------------------------------------!
107	! FLOAT INITILISATION: NO RESTART FILE !
108	!--------------------------------------!
109	ELSE !ln_rstflo
110
111	IF( ln_ariane )THEN !Add new floats with ariane convention
112	CALL flo_add_new_ariane_floats(1,jpnfl)
113	ELSE !Add new floats with long/lat convention
114	CALL flo_add_new_floats(1,jpnfl)
115	ENDIF
116
117	ENDIF
118
119	END SUBROUTINE flo_dom
120
121	SUBROUTINE flo_add_new_floats(kfl_start, kfl_end)
122	!! -------------------------------------------------------------
123	!! * SUBROUTINE add_new_arianefloats *
124	!!
125	!! ** Purpose :
126	!!
127	!! First initialisation of floats
128	!! the initials positions of floats are written in a file
129	!! with a variable to know if it is a isobar float a number
130	!! to identified who want the trajectories of this float and
131	!! an index for the number of the float
132	!! open the init file
133	!!
134	!! ** Method :
135	!!----------------------------------------------------------------------
136	INTEGER, INTENT(in) :: kfl_start, kfl_end
137	!!
138	INTEGER :: inum ! file unit
139	INTEGER :: jfl,ji, jj, jk ! dummy loop indices
140	INTEGER :: itrash ! trash var for reading
141	INTEGER :: ifl ! number of floats to read
142	REAL(wp) :: zdxab, zdyad
143	LOGICAL :: llinmesh
144	CHARACTER(len=80) :: cltmp
145	!!---------------------------------------------------------------------
146	ifl = kfl_end-kfl_start+1
147
148	! we get the init values
149	!-----------------------
150	CALL ctl_opn( inum , clname1, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
151	DO jfl = kfl_start,kfl_end
152	READ(inum,*) flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),itrash
153	IF(lwp) THEN
154	write(numout,*)'read:',jfl,flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),itrash
155	IF(lflush) CALL FLUSH(numout)
156	ENDIF
157	END DO
158	CLOSE(inum)
159
160	! Test to find the grid point coordonate with the geographical position
161	!----------------------------------------------------------------------
162	DO jfl = kfl_start,kfl_end
163	ihtest(jfl) = 0
164	ivtest(jfl) = 0
165	ikmfl(jfl) = 0
166	# if defined key_mpp_mpi
167	DO ji = MAX(nldi,2), nlei
168	DO jj = MAX(nldj,2), nlej ! NO vector opt.
169	# else
170	DO ji = 2, jpi
171	DO jj = 2, jpj ! NO vector opt.
172	# endif
173	! For each float we find the indexes of the mesh
174	CALL flo_findmesh(glamf(ji-1,jj-1),gphif(ji-1,jj-1), &
175	glamf(ji-1,jj ),gphif(ji-1,jj ), &
176	glamf(ji ,jj ),gphif(ji ,jj ), &
177	glamf(ji ,jj-1),gphif(ji ,jj-1), &
178	flxx(jfl) ,flyy(jfl) , &
179	glamt(ji ,jj ),gphit(ji ,jj ), llinmesh)
180	IF( llinmesh )THEN
181	iimfl(jfl) = ji
182	ijmfl(jfl) = jj
183	ihtest(jfl) = ihtest(jfl)+1
184	DO jk = 1, jpk-1
185	IF( (gdepw_n(ji,jj,jk) <= flzz(jfl)) .AND. (gdepw_n(ji,jj,jk+1) > flzz(jfl)) ) THEN
186	ikmfl(jfl) = jk
187	ivtest(jfl) = ivtest(jfl) + 1
188	ENDIF
189	END DO
190	ENDIF
191	END DO
192	END DO
193
194	! If the float is in a mesh computed by an other processor we put iimfl=ijmfl=-1
195	IF( ihtest(jfl) == 0 ) THEN
196	iimfl(jfl) = -1
197	ijmfl(jfl) = -1
198	ENDIF
199	END DO
200
201	!Test if each float is in one and only one proc
202	!----------------------------------------------
203	IF( lk_mpp ) THEN
204	CALL mpp_sum('flodom', ihtest,jpnfl)
205	CALL mpp_sum('flodom', ivtest,jpnfl)
206	ENDIF
207	DO jfl = kfl_start,kfl_end
208
209	IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1) ) THEN
210	WRITE(cltmp,'(A10,i4.4,A20)' )'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH'
211	CALL ctl_stop('STOP',TRIM(cltmp) )
212	ENDIF
213	IF( (ihtest(jfl) == 0) ) THEN
214	WRITE(cltmp,'(A10,i4.4,A20)' )'THE FLOAT',jfl,' IS IN NO MESH'
215	CALL ctl_stop('STOP',TRIM(cltmp) )
216	ENDIF
217	END DO
218
219	! We compute the distance between the float and the face of the mesh
220	!-------------------------------------------------------------------
221	DO jfl = kfl_start,kfl_end
222
223	! Made only if the float is in the domain of the processor
224	IF( (iimfl(jfl) >= 0) .AND. (ijmfl(jfl) >= 0) ) THEN
225
226	! TEST TO KNOW IF THE FLOAT IS NOT INITIALISED IN THE COAST
227	idomfl(jfl) = 0
228	IF( tmask(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) == 0. ) idomfl(jfl) = 1
229
230	! Computation of the distance between the float and the faces of the mesh
231	! zdxab
232	! .
233	! B----.---------C
234	! \| . \|
235	! \|<------>flo \|
236	! \| ^ \|
237	! \| \|.....\|....zdyad
238	! \| \| \|
239	! A--------\|-----D
240	!
241	zdxab = flo_dstnce( flxx(jfl), flyy(jfl), glamf(iimfl(jfl)-1,ijmfl(jfl)-1), flyy(jfl) )
242	zdyad = flo_dstnce( flxx(jfl), flyy(jfl), flxx(jfl), gphif(iimfl(jfl)-1,ijmfl(jfl)-1) )
243
244	! Translation of this distances (in meter) in indexes
245	zgifl(jfl)= (iimfl(jfl)-0.5) + zdxab/e1u(iimfl(jfl)-1,ijmfl(jfl)) + (mig(1)-1)
246	zgjfl(jfl)= (ijmfl(jfl)-0.5) + zdyad/e2v(iimfl(jfl),ijmfl(jfl)-1) + (mjg(1)-1)
247	zgkfl(jfl) = (( gdepw_n(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) - flzz(jfl) )* ikmfl(jfl)) &
248	& / ( gdepw_n(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) &
249	& - gdepw_n(iimfl(jfl),ijmfl(jfl),ikmfl(jfl) ) ) &
250	& + (( flzz(jfl)-gdepw_n(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) ) *(ikmfl(jfl)+1)) &
251	& / ( gdepw_n(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)+1) &
252	& - gdepw_n(iimfl(jfl),ijmfl(jfl),ikmfl(jfl)) )
253	ELSE
254	zgifl(jfl) = 0.e0
255	zgjfl(jfl) = 0.e0
256	zgkfl(jfl) = 0.e0
257	ENDIF
258
259	END DO
260
261	! The sum of all the arrays zgifl, zgjfl, zgkfl give 3 arrays with the positions of all the floats.
262	IF( lk_mpp ) THEN
263	CALL mpp_sum( 'flodom', zgjfl, ifl ) ! sums over the global domain
264	CALL mpp_sum( 'flodom', zgkfl, ifl )
265	ENDIF
266
267	DO jfl = kfl_start,kfl_end
268	tpifl(jfl) = zgifl(jfl)
269	tpjfl(jfl) = zgjfl(jfl)
270	tpkfl(jfl) = zgkfl(jfl)
271	END DO
272
273	! WARNING : initial position not in the sea
274	IF( .NOT. ln_rstflo ) THEN
275	DO jfl = kfl_start,kfl_end
276	IF( idomfl(jfl) == 1 ) THEN
277	IF(lwp) THEN
278	WRITE(numout,)'****************************'
279	WRITE(numout,*)'!!!!!!! WARNING !!!!!!!!!!'
280	WRITE(numout,)'****************************'
281	WRITE(numout,*)'The float number',jfl,'is out of the sea.'
282	WRITE(numout,*)'geographical position',flxx(jfl),flyy(jfl),flzz(jfl)
283	WRITE(numout,*)'index position',tpifl(jfl),tpjfl(jfl),tpkfl(jfl)
284	IF(lflush) CALL FLUSH(numout)
285	ENDIF
286	ENDIF
287	END DO
288	ENDIF
289
290	END SUBROUTINE flo_add_new_floats
291
292	SUBROUTINE flo_add_new_ariane_floats(kfl_start, kfl_end)
293	!! -------------------------------------------------------------
294	!! * SUBROUTINE add_new_arianefloats *
295	!!
296	!! ** Purpose :
297	!! First initialisation of floats with ariane convention
298	!!
299	!! The indexes are read directly from file (warning ariane
300	!! convention, are refered to
301	!! U,V,W grids - and not T-)
302	!! The isobar advection is managed with the sign of tpkfl ( >0 -> 3D
303	!! advection, <0 -> 2D)
304	!! Some variables are not read, as - gl : time index; 4th
305	!! column
306	!! - transport : transport ; 5th
307	!! column
308	!! and paste in the jtrash var
309	!! At the end, ones need to replace the indexes on T grid
310	!! RMQ : there is no float groups identification !
311	!!
312	!!
313	!! ** Method :
314	!!----------------------------------------------------------------------
315	INTEGER, INTENT(in) :: kfl_start, kfl_end
316	!!
317	INTEGER :: inum ! file unit
318	INTEGER :: ierr, ifl
319	INTEGER :: jfl, jfl1 ! dummy loop indices
320	INTEGER :: itrash ! trash var for reading
321	CHARACTER(len=80) :: cltmp
322
323	!!----------------------------------------------------------------------
324	nisobfl(kfl_start:kfl_end) = 1 ! we assume that by default we want 3D advection
325
326	ifl = kfl_end - kfl_start + 1 ! number of floats to read
327
328	! we check that the number of floats in the init_file are consistant with the namelist
329	IF( lwp ) THEN
330
331	jfl1=0
332	ierr=0
333	CALL ctl_opn( inum, clname2, 'OLD', 'FORMATTED', 'SEQUENTIAL', 1, numout, .TRUE., 1 )
334	DO WHILE (ierr .EQ. 0)
335	jfl1=jfl1+1
336	READ(inum,*, iostat=ierr)
337	END DO
338	CLOSE(inum)
339	IF( (jfl1-1) .NE. ifl )THEN
340	WRITE(cltmp,'(A25,A20,A3,i4.4,A10,i4.4)')"the number of floats in ",TRIM(clname2), &
341	" = ",jfl1," is not equal to jfl= ",ifl
342	CALL ctl_stop('STOP',TRIM(cltmp) )
343	ENDIF
344
345	ENDIF
346
347	! we get the init values
348	CALL ctl_opn( inum, clname2, 'OLD', 'FORMATTED', 'SEQUENTIAL', 1, numout, .TRUE., 1 )
349	DO jfl = kfl_start, kfl_end
350	READ(inum,*) tpifl(jfl),tpjfl(jfl),tpkfl(jfl),itrash, itrash
351
352	IF ( tpkfl(jfl) .LT. 0. ) nisobfl(jfl) = 0 !set the 2D advection according to init_float
353	ngrpfl(jfl)=jfl
354	END DO
355
356	! conversion from ariane index to T grid index
357	tpkfl(kfl_start:kfl_end) = abs(tpkfl)-0.5 ! reversed vertical axis
358	tpifl(kfl_start:kfl_end) = tpifl+0.5
359	tpjfl(kfl_start:kfl_end) = tpjfl+0.5
360
361
362	END SUBROUTINE flo_add_new_ariane_floats
363
364
365	SUBROUTINE flo_findmesh( pax, pay, pbx, pby, &
366	pcx, pcy, pdx, pdy, &
367	px ,py ,ptx, pty, ldinmesh )
368	!! -------------------------------------------------------------
369	!! * ROUTINE findmesh *
370	!!
371	!! ** Purpose : Find the index of mesh for the point spx spy
372	!!
373	!! ** Method :
374	!!----------------------------------------------------------------------
375	REAL(wp) :: &
376	pax, pay, pbx, pby, & ! ???
377	pcx, pcy, pdx, pdy, & ! ???
378	px, py, & ! longitude and latitude
379	ptx, pty ! ???
380	LOGICAL :: ldinmesh ! ???
381	!!
382	REAL(wp) :: zabt, zbct, zcdt, zdat, zabpt, zbcpt, zcdpt, zdapt
383	!!---------------------------------------------------------------------
384	!! Statement function
385	REAL(wp) :: fsline
386	REAL(wp) :: psax, psay, psbx, psby, psx, psy
387	fsline( psax, psay, psbx, psby, psx, psy ) = psy * ( psbx - psax ) &
388	& - psx * ( psby - psay ) &
389	& + psax * psby - psay * psbx
390	!!---------------------------------------------------------------------
391
392	! 4 semi plane defined by the 4 points and including the T point
393	zabt = fsline(pax,pay,pbx,pby,ptx,pty)
394	zbct = fsline(pbx,pby,pcx,pcy,ptx,pty)
395	zcdt = fsline(pcx,pcy,pdx,pdy,ptx,pty)
396	zdat = fsline(pdx,pdy,pax,pay,ptx,pty)
397
398	! 4 semi plane defined by the 4 points and including the extrememity
399	zabpt = fsline(pax,pay,pbx,pby,px,py)
400	zbcpt = fsline(pbx,pby,pcx,pcy,px,py)
401	zcdpt = fsline(pcx,pcy,pdx,pdy,px,py)
402	zdapt = fsline(pdx,pdy,pax,pay,px,py)
403
404	! We compare the semi plane T with the semi plane including the point
405	! to know if it is in this mesh.
406	! For numerical reasons it is possible that for a point which is on
407	! the line we don't have exactly zero with fsline function. We want
408	! that a point can't be in 2 mesh in the same time, so we put the
409	! coefficient to zero if it is smaller than 1.E-12
410
411	IF( ABS(zabpt) <= 1.E-12 ) zabpt = 0.
412	IF( ABS(zbcpt) <= 1.E-12 ) zbcpt = 0.
413	IF( ABS(zcdpt) <= 1.E-12 ) zcdpt = 0.
414	IF( ABS(zdapt) <= 1.E-12 ) zdapt = 0.
415	IF( (zabtzabpt > 0.) .AND. (zbctzbcpt >= 0. ) .AND. ( zcdtzcdpt >= 0. ) .AND. ( zdatzdapt > 0. ) &
416	.AND. ( px <= MAX(pcx,pdx) ) .AND. ( px >= MIN(pax,pbx) ) &
417	.AND. ( py <= MAX(pby,pcy) ) .AND. ( py >= MIN(pay,pdy) ) ) THEN
418	ldinmesh=.TRUE.
419	ELSE
420	ldinmesh=.FALSE.
421	ENDIF
422	!
423	END SUBROUTINE flo_findmesh
424
425
426	FUNCTION flo_dstnce( pla1, phi1, pla2, phi2 )
427	!! -------------------------------------------------------------
428	!! * Function dstnce *
429	!!
430	!! ** Purpose : returns distance (in m) between two geographical
431	!! points
432	!! ** Method :
433	!!----------------------------------------------------------------------
434	REAL(wp), INTENT(in) :: pla1, phi1, pla2, phi2 ! ???
435	!!
436	REAL(wp) :: dly1, dly2, dlx1, dlx2, dlx, dls, dld, dpi
437	REAL(wp) :: flo_dstnce
438	!!---------------------------------------------------------------------
439	!
440	dpi = 2._wp * ASIN(1._wp)
441	dls = dpi / 180._wp
442	dly1 = phi1 * dls
443	dly2 = phi2 * dls
444	dlx1 = pla1 * dls
445	dlx2 = pla2 * dls
446	!
447	dlx = SIN(dly1) * SIN(dly2) + COS(dly1) * COS(dly2) * COS(dlx2-dlx1)
448	!
449	IF( ABS(dlx) > 1.0_wp ) dlx = 1.0_wp
450	!
451	dld = ATAN(DSQRT( 1._wp * ( 1._wp-dlx )/( 1._wp+dlx ) )) * 222.24_wp / dls
452	flo_dstnce = dld * 1000._wp
453	!
454	END FUNCTION flo_dstnce
455
456	INTEGER FUNCTION flo_dom_alloc()
457	!!----------------------------------------------------------------------
458	!! * FUNCTION flo_dom_alloc *
459	!!----------------------------------------------------------------------
460
461	ALLOCATE( iimfl(jpnfl) , ijmfl(jpnfl) , ikmfl(jpnfl) , &
462	idomfl(jpnfl), ivtest(jpnfl), ihtest(jpnfl), &
463	zgifl(jpnfl) , zgjfl(jpnfl) , zgkfl(jpnfl) , STAT=flo_dom_alloc )
464	!
465	CALL mpp_sum ( 'flodom', flo_dom_alloc )
466	IF( flo_dom_alloc /= 0 ) CALL ctl_stop( 'STOP', 'flo_dom_alloc: failed to allocate arrays' )
467	END FUNCTION flo_dom_alloc
468
469
470	#else
471	!!----------------------------------------------------------------------
472	!! Default option Empty module
473	!!----------------------------------------------------------------------
474	CONTAINS
475	SUBROUTINE flo_dom ! Empty routine
476	WRITE(,) 'flo_dom: : You should not have seen this print! error?'
477	END SUBROUTINE flo_dom
478	#endif
479
480	!!======================================================================
481	END MODULE flodom

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/FLO/flodom.F90 @ 10986

Download in other formats: