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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

flodom.F90 in NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/FLO – NEMO

Context Navigation

source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/FLO/flodom.F90 @ 12165

Last change on this file since 12165 was 12165, checked in by davestorkey, 4 years ago
2019/dev_ASINTER-01-05_merged: Update to r12072 of trunk.
Property svn:keywords set to `Id`
File size: 19.5 KB

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

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

Download in other formats: