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.

icbutl.F90 in NEMO/branches/2019/fix_ticket2238_solution1/src/OCE/ICB – NEMO

Context Navigation

source: NEMO/branches/2019/fix_ticket2238_solution1/src/OCE/ICB/icbutl.F90 @ 10677

Last change on this file since 10677 was 10677, checked in by mathiot, 5 years ago
issue for EW reproducibility
Property svn:keywords set to `Id`
File size: 34.3 KB

Line
1	MODULE icbutl
2	!!======================================================================
3	!! * MODULE icbutl *
4	!! Icebergs: various iceberg utility routines
5	!!======================================================================
6	!! History : 3.3.1 ! 2010-01 (Martin&Adcroft) Original code
7	!! - ! 2011-03 (Madec) Part conversion to NEMO form
8	!! - ! Removal of mapping from another grid
9	!! - ! 2011-04 (Alderson) Split into separate modules
10	!!----------------------------------------------------------------------
11
12	!!----------------------------------------------------------------------
13	!! icb_utl_interp :
14	!! icb_utl_bilin :
15	!! icb_utl_bilin_e :
16	!!----------------------------------------------------------------------
17	USE par_oce ! ocean parameters
18	USE dom_oce ! ocean domain
19	USE in_out_manager ! IO parameters
20	USE lbclnk ! lateral boundary condition
21	USE lib_mpp ! MPI code and lk_mpp in particular
22	USE icb_oce ! define iceberg arrays
23	USE sbc_oce ! ocean surface boundary conditions
24	#if defined key_si3
25	USE ice, ONLY: u_ice, v_ice, hm_i ! SI3 variables
26	USE icevar ! ice_var_sshdyn
27	USE sbc_ice, ONLY: snwice_mass, snwice_mass_b
28	#endif
29
30	IMPLICIT NONE
31	PRIVATE
32
33	PUBLIC icb_utl_copy ! routine called in icbstp module
34	PUBLIC icb_utl_interp ! routine called in icbdyn, icbthm modules
35	PUBLIC icb_utl_bilin ! routine called in icbini, icbdyn modules
36	PUBLIC icb_utl_bilin_x ! routine called in icbdyn module
37	PUBLIC icb_utl_add ! routine called in icbini.F90, icbclv, icblbc and icbrst modules
38	PUBLIC icb_utl_delete ! routine called in icblbc, icbthm modules
39	PUBLIC icb_utl_destroy ! routine called in icbstp module
40	PUBLIC icb_utl_track ! routine not currently used, retain just in case
41	PUBLIC icb_utl_print_berg ! routine called in icbthm module
42	PUBLIC icb_utl_print ! routine called in icbini, icbstp module
43	PUBLIC icb_utl_count ! routine called in icbdia, icbini, icblbc, icbrst modules
44	PUBLIC icb_utl_incr ! routine called in icbini, icbclv modules
45	PUBLIC icb_utl_yearday ! routine called in icbclv, icbstp module
46	PUBLIC icb_utl_mass ! routine called in icbdia module
47	PUBLIC icb_utl_heat ! routine called in icbdia module
48
49	!!----------------------------------------------------------------------
50	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
51	!! $Id$
52	!! Software governed by the CeCILL license (see ./LICENSE)
53	!!----------------------------------------------------------------------
54	CONTAINS
55
56	SUBROUTINE icb_utl_copy()
57	!!----------------------------------------------------------------------
58	!! * ROUTINE icb_utl_copy *
59	!!
60	!! ** Purpose : iceberg initialization.
61	!!
62	!! ** Method : - blah blah
63	!!----------------------------------------------------------------------
64	#if defined key_si3
65	REAL(wp), DIMENSION(jpi,jpj) :: zssh_lead_m ! ocean surface (ssh_m) if ice is not embedded
66	! ! ocean surface in leads if ice is embedded
67	#endif
68	! copy nemo forcing arrays into iceberg versions with extra halo
69	! only necessary for variables not on T points
70	! and ssh which is used to calculate gradients
71
72	uo_e(:,:) = 0._wp ; uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1)
73	vo_e(:,:) = 0._wp ; vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1)
74	ff_e(:,:) = 0._wp ; ff_e(1:jpi,1:jpj) = ff_f (:,:)
75	tt_e(:,:) = 0._wp ; tt_e(1:jpi,1:jpj) = sst_m(:,:)
76	fr_e(:,:) = 0._wp ; fr_e(1:jpi,1:jpj) = fr_i (:,:)
77	ua_e(:,:) = 0._wp ; ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
78	va_e(:,:) = 0._wp ; va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk
79	!
80	CALL lbc_lnk_icb( 'icbutl', uo_e, 'U', -1._wp, 1, 1 )
81	CALL lbc_lnk_icb( 'icbutl', vo_e, 'V', -1._wp, 1, 1 )
82	CALL lbc_lnk_icb( 'icbutl', ff_e, 'F', +1._wp, 1, 1 )
83	CALL lbc_lnk_icb( 'icbutl', ua_e, 'U', -1._wp, 1, 1 )
84	CALL lbc_lnk_icb( 'icbutl', va_e, 'V', -1._wp, 1, 1 )
85	CALL lbc_lnk_icb( 'icbutl', fr_e, 'T', +1._wp, 1, 1 )
86	CALL lbc_lnk_icb( 'icbutl', tt_e, 'T', +1._wp, 1, 1 )
87	#if defined key_si3
88	hicth(:,:) = 0._wp ; hicth(1:jpi,1:jpj) = hm_i (:,:)
89	ui_e(:,:) = 0._wp ; ui_e(1:jpi, 1:jpj) = u_ice(:,:)
90	vi_e(:,:) = 0._wp ; vi_e(1:jpi, 1:jpj) = v_ice(:,:)
91	!
92	! compute ssh slope using ssh_lead if embedded
93	zssh_lead_m(:,:) = ice_var_sshdyn(ssh_m, snwice_mass, snwice_mass_b)
94	ssh_e(:,:) = 0._wp ; ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1)
95	!
96	CALL lbc_lnk_icb( 'icbutl', hicth, 'T', +1._wp, 1, 1 )
97	CALL lbc_lnk_icb( 'icbutl', ui_e , 'U', -1._wp, 1, 1 )
98	CALL lbc_lnk_icb( 'icbutl', vi_e , 'V', -1._wp, 1, 1 )
99	#else
100	ssh_e(:,:) = 0._wp ; ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)
101	#endif
102	CALL lbc_lnk_icb( 'icbutl', ssh_e, 'T', +1._wp, 1, 1 )
103	!
104	END SUBROUTINE icb_utl_copy
105
106
107	SUBROUTINE icb_utl_interp( pi, pe1, puo, pui, pua, pssh_i, &
108	& pj, pe2, pvo, pvi, pva, pssh_j, &
109	& psst, pcn, phi, pff )
110	!!----------------------------------------------------------------------
111	!! * ROUTINE icb_utl_interp *
112	!!
113	!! ** Purpose : interpolation
114	!!
115	!! ** Method : - interpolate from various ocean arrays onto iceberg position
116	!!
117	!! !!gm CAUTION here I do not care of the slip/no-slip conditions
118	!! this can be done later (not that easy to do...)
119	!! right now, U is 0 in land so that the coastal value of velocity parallel to the coast
120	!! is half the off shore value, wile the normal-to-the-coast value is zero.
121	!! This is OK as a starting point.
122	!!
123	!!----------------------------------------------------------------------
124	REAL(wp), INTENT(in ) :: pi , pj ! position in (i,j) referential
125	REAL(wp), INTENT( out) :: pe1, pe2 ! i- and j scale factors
126	REAL(wp), INTENT( out) :: puo, pvo, pui, pvi, pua, pva ! ocean, ice and wind speeds
127	REAL(wp), INTENT( out) :: pssh_i, pssh_j ! ssh i- & j-gradients
128	REAL(wp), INTENT( out) :: psst, pcn, phi, pff ! SST, ice concentration, ice thickness, Coriolis
129	!
130	REAL(wp) :: zcd, zmod ! local scalars
131	!!----------------------------------------------------------------------
132
133	pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) ! scale factors
134	pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj )
135	!
136	puo = icb_utl_bilin_h( uo_e, pi, pj, 'U' ) ! ocean velocities
137	pvo = icb_utl_bilin_h( vo_e, pi, pj, 'V' )
138	psst = icb_utl_bilin_h( tt_e, pi, pj, 'T' ) ! SST
139	pcn = icb_utl_bilin_h( fr_e , pi, pj, 'T' ) ! ice concentration
140	pff = icb_utl_bilin_h( ff_e , pi, pj, 'F' ) ! Coriolis parameter
141	!
142	pua = icb_utl_bilin_h( ua_e , pi, pj, 'U' ) ! 10m wind
143	pva = icb_utl_bilin_h( va_e , pi, pj, 'V' ) ! here (ua,va) are stress => rough conversion from stress to speed
144	zcd = 1.22_wp * 1.5e-3_wp ! air density * drag coefficient
145	zmod = 1._wp / MAX( 1.e-20, SQRT( zcd * SQRT( puapua + pvapva) ) )
146	pua = pua * zmod ! note: stress module=0 necessarly implies ua=va=0
147	pva = pva * zmod
148
149	#if defined key_si3
150	pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' ) ! sea-ice velocities
151	pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' )
152	phi = icb_utl_bilin_h( hicth, pi, pj, 'T' ) ! ice thickness
153	#else
154	pui = 0._wp
155	pvi = 0._wp
156	phi = 0._wp
157	#endif
158
159	! Estimate SSH gradient in i- and j-direction (centred evaluation)
160	pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T' ) - &
161	& icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' ) ) / ( 0.2_wp * pe1 )
162	pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T' ) - &
163	& icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' ) ) / ( 0.2_wp * pe2 )
164	!
165	END SUBROUTINE icb_utl_interp
166
167
168	REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type )
169	!!----------------------------------------------------------------------
170	!! * FUNCTION icb_utl_bilin *
171	!!
172	!! ** Purpose : bilinear interpolation at berg location depending on the grid-point type
173	!! this version deals with extra halo points
174	!!
175	!! !!gm CAUTION an optional argument should be added to handle
176	!! the slip/no-slip conditions ==>>> to be done later
177	!!
178	!!----------------------------------------------------------------------
179	REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) :: pfld ! field to be interpolated
180	REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential
181	CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points
182	!
183	INTEGER :: ii, ij ! local integer
184	REAL(wp) :: zi, zj ! local real
185	!!----------------------------------------------------------------------
186	!
187	SELECT CASE ( cd_type )
188	CASE ( 'T' )
189	! note that here there is no +0.5 added
190	! since we're looking for four T points containing quadrant we're in of
191	! current T cell
192	ii = MAX(0, INT( pi ))
193	ij = MAX(0, INT( pj )) ! T-point
194	zi = pi - REAL(ii,wp)
195	zj = pj - REAL(ij,wp)
196	CASE ( 'U' )
197	ii = MAX(0, INT( pi-0.5_wp ))
198	ij = MAX(0, INT( pj )) ! U-point
199	zi = pi - 0.5_wp - REAL(ii,wp)
200	zj = pj - REAL(ij,wp)
201	CASE ( 'V' )
202	ii = MAX(0, INT( pi ))
203	ij = MAX(0, INT( pj-0.5_wp )) ! V-point
204	zi = pi - REAL(ii,wp)
205	zj = pj - 0.5_wp - REAL(ij,wp)
206	CASE ( 'F' )
207	ii = MAX(0, INT( pi-0.5_wp ))
208	ij = MAX(0, INT( pj-0.5_wp )) ! F-point
209	zi = pi - 0.5_wp - REAL(ii,wp)
210	zj = pj - 0.5_wp - REAL(ij,wp)
211	END SELECT
212	!
213	! find position in this processor. Prevent near edge problems (see #1389)
214	! (PM) will be useless if extra halo is used in NEMO
215	!
216	IF ( ii <= mig(1)-1 ) THEN ; ii = 0
217	ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi
218	ELSE ; ii = mi1(ii)
219	ENDIF
220	IF ( ij <= mjg(1)-1 ) THEN ; ij = 0
221	ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj
222	ELSE ; ij = mj1(ij)
223	ENDIF
224	!
225	!
226	icb_utl_bilin_h = ( pfld(ii,ij ) * (1._wp-zi) + pfld(ii+1,ij ) * zi ) * (1._wp-zj) &
227	& + ( pfld(ii,ij+1) * (1._wp-zi) + pfld(ii+1,ij+1) * zi ) * zj
228	!
229	END FUNCTION icb_utl_bilin_h
230
231
232	REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type )
233	!!----------------------------------------------------------------------
234	!! * FUNCTION icb_utl_bilin *
235	!!
236	!! ** Purpose : bilinear interpolation at berg location depending on the grid-point type
237	!!
238	!! !!gm CAUTION an optional argument should be added to handle
239	!! the slip/no-slip conditions ==>>> to be done later
240	!!
241	!!----------------------------------------------------------------------
242	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated
243	REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential
244	CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points
245	!
246	INTEGER :: ii, ij ! local integer
247	REAL(wp) :: zi, zj ! local real
248	!!----------------------------------------------------------------------
249	!
250	SELECT CASE ( cd_type )
251	CASE ( 'T' )
252	! note that here there is no +0.5 added
253	! since we're looking for four T points containing quadrant we're in of
254	! current T cell
255	ii = MAX(1, INT( pi ))
256	ij = MAX(1, INT( pj )) ! T-point
257	zi = pi - REAL(ii,wp)
258	zj = pj - REAL(ij,wp)
259	CASE ( 'U' )
260	ii = MAX(1, INT( pi-0.5 ))
261	ij = MAX(1, INT( pj )) ! U-point
262	zi = pi - 0.5 - REAL(ii,wp)
263	zj = pj - REAL(ij,wp)
264	CASE ( 'V' )
265	ii = MAX(1, INT( pi ))
266	ij = MAX(1, INT( pj-0.5 )) ! V-point
267	zi = pi - REAL(ii,wp)
268	zj = pj - 0.5 - REAL(ij,wp)
269	CASE ( 'F' )
270	ii = MAX(1, INT( pi-0.5 ))
271	ij = MAX(1, INT( pj-0.5 )) ! F-point
272	zi = pi - 0.5 - REAL(ii,wp)
273	zj = pj - 0.5 - REAL(ij,wp)
274	END SELECT
275	!
276	! find position in this processor. Prevent near edge problems (see #1389)
277	IF ( ii < mig( 1 ) ) THEN ; ii = 1
278	ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi
279	ELSE ; ii = mi1(ii)
280	ENDIF
281	IF ( ij < mjg( 1 ) ) THEN ; ij = 1
282	ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj
283	ELSE ; ij = mj1(ij)
284	ENDIF
285	!
286	IF( ii == jpi ) ii = ii-1
287	IF( ij == jpj ) ij = ij-1
288	!
289	icb_utl_bilin = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) &
290	& + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj
291	!
292	END FUNCTION icb_utl_bilin
293
294
295	REAL(wp) FUNCTION icb_utl_bilin_x( pfld, pi, pj )
296	!!----------------------------------------------------------------------
297	!! * FUNCTION icb_utl_bilin_x *
298	!!
299	!! ** Purpose : bilinear interpolation at berg location depending on the grid-point type
300	!! Special case for interpolating longitude
301	!!
302	!! !!gm CAUTION an optional argument should be added to handle
303	!! the slip/no-slip conditions ==>>> to be done later
304	!!
305	!!----------------------------------------------------------------------
306	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated
307	REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential
308	!
309	INTEGER :: ii, ij ! local integer
310	REAL(wp) :: zi, zj ! local real
311	REAL(wp) :: zret ! local real
312	REAL(wp), DIMENSION(4) :: z4
313	!!----------------------------------------------------------------------
314	!
315	! note that here there is no +0.5 added
316	! since we're looking for four T points containing quadrant we're in of
317	! current T cell
318	ii = MAX(1, INT( pi ))
319	ij = MAX(1, INT( pj )) ! T-point
320	zi = pi - REAL(ii,wp)
321	zj = pj - REAL(ij,wp)
322	!
323	! find position in this processor. Prevent near edge problems (see #1389)
324	IF ( ii < mig( 1 ) ) THEN ; ii = 1
325	ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi
326	ELSE ; ii = mi1(ii)
327	ENDIF
328	IF ( ij < mjg( 1 ) ) THEN ; ij = 1
329	ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj
330	ELSE ; ij = mj1(ij)
331	ENDIF
332	!
333	IF( ii == jpi ) ii = ii-1
334	IF( ij == jpj ) ij = ij-1
335	!
336	z4(1) = pfld(ii ,ij )
337	z4(2) = pfld(ii+1,ij )
338	z4(3) = pfld(ii ,ij+1)
339	z4(4) = pfld(ii+1,ij+1)
340	IF( MAXVAL(z4) - MINVAL(z4) > 90._wp ) THEN
341	WHERE( z4 < 0._wp ) z4 = z4 + 360._wp
342	ENDIF
343	!
344	zret = (z4(1) * (1.-zi) + z4(2) * zi) * (1.-zj) + (z4(3) * (1.-zi) + z4(4) * zi) * zj
345	IF( zret > 180._wp ) zret = zret - 360._wp
346	icb_utl_bilin_x = zret
347	!
348	END FUNCTION icb_utl_bilin_x
349
350
351	REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj )
352	!!----------------------------------------------------------------------
353	!! * FUNCTION dom_init *
354	!!
355	!! ** Purpose : bilinear interpolation at berg location of horizontal scale factor
356	!! ** Method : interpolation done using the 4 nearest grid points among
357	!! t-, u-, v-, and f-points.
358	!!----------------------------------------------------------------------
359	REAL(wp), DIMENSION(:,:), INTENT(in) :: pet, peu, pev, pef ! horizontal scale factor to be interpolated at t-,u-,v- & f-pts
360	REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential
361	!
362	INTEGER :: ii, ij, icase ! local integer
363	!
364	! weights corresponding to corner points of a T cell quadrant
365	REAL(wp) :: zi, zj ! local real
366	!
367	! values at corner points of a T cell quadrant
368	! 00 = bottom left, 10 = bottom right, 01 = top left, 11 = top right
369	REAL(wp) :: ze00, ze10, ze01, ze11
370	!!----------------------------------------------------------------------
371	!
372	ii = MAX(1, INT( pi )) ; ij = MAX(1, INT( pj )) ! left bottom T-point (i,j) indices
373
374	! fractional box spacing
375	! 0 <= zi < 0.5, 0 <= zj < 0.5 --> NW quadrant of current T cell
376	! 0.5 <= zi < 1 , 0 <= zj < 0.5 --> NE quadrant
377	! 0 <= zi < 0.5, 0.5 <= zj < 1 --> SE quadrant
378	! 0.5 <= zi < 1 , 0.5 <= zj < 1 --> SW quadrant
379
380	zi = pi - REAL(ii,wp) !!gm use here mig, mjg arrays
381	zj = pj - REAL(ij,wp)
382
383	! find position in this processor. Prevent near edge problems (see #1389)
384	IF ( ii < mig( 1 ) ) THEN ; ii = 1 ; PRINT *, 'SHOULD NOT SEE THIS 1'
385	ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi ; PRINT *, 'SHOULD NOT SEE THIS jpi'
386	ELSE ; ii = mi1(ii)
387	ENDIF
388	IF ( ij < mjg( 1 ) ) THEN ; ij = 1 ; PRINT *, 'SHOULD NOT SEE THIS i 1'
389	ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj ; PRINT *, 'SHOULD NOT SEE THIS jpj'
390	ELSE ; ij = mj1(ij)
391	ENDIF
392	!
393	IF( ii == jpi ) THEN ; ii = ii-1 ; PRINT *, 'SHOULD NOT SEE THIS ii-1' ; END IF
394	IF( ij == jpj ) THEN ; ij = ij-1 ; PRINT *, 'SHOULD NOT SEE THIS ij-1' ; END IF
395	!
396	IF( 0.0_wp <= zi .AND. zi < 0.5_wp ) THEN
397	IF( 0.0_wp <= zj .AND. zj < 0.5_wp ) THEN ! NE quadrant
398	! ! i=I i=I+1/2
399	ze01 = pev(ii ,ij ) ; ze11 = pef(ii ,ij ) ! j=J+1/2 V ------- F
400	ze00 = pet(ii ,ij ) ; ze10 = peu(ii ,ij ) ! j=J T ------- U
401	zi = 2._wp * zi
402	zj = 2._wp * zj
403	ELSE ! SE quadrant
404	! ! i=I i=I+1/2
405	ze01 = pet(ii ,ij+1) ; ze11 = peu(ii ,ij+1) ! j=J+1 T ------- U
406	ze00 = pev(ii ,ij ) ; ze10 = pef(ii ,ij ) ! j=J+1/2 V ------- F
407	zi = 2._wp * zi
408	zj = 2._wp * (zj-0.5_wp)
409	ENDIF
410	ELSE
411	IF( 0.0_wp <= zj .AND. zj < 0.5_wp ) THEN ! NW quadrant
412	! ! i=I i=I+1/2
413	ze01 = pef(ii ,ij ) ; ze11 = pev(ii+1,ij) ! j=J+1/2 F ------- V
414	ze00 = peu(ii ,ij ) ; ze10 = pet(ii+1,ij) ! j=J U ------- T
415	zi = 2._wp * (zi-0.5_wp)
416	zj = 2._wp * zj
417	ELSE ! SW quadrant
418	! ! i=I+1/2 i=I+1
419	ze01 = peu(ii ,ij+1) ; ze11 = pet(ii+1,ij+1) ! j=J+1 U ------- T
420	ze00 = pef(ii ,ij ) ; ze10 = pev(ii+1,ij ) ! j=J+1/2 F ------- V
421	zi = 2._wp * (zi-0.5_wp)
422	zj = 2._wp * (zj-0.5_wp)
423	ENDIF
424	ENDIF
425	!
426	icb_utl_bilin_e = ( ze01 * (1._wp-zi) + ze11 * zi ) * zj &
427	& + ( ze00 * (1._wp-zi) + ze10 * zi ) * (1._wp-zj)
428	!
429	END FUNCTION icb_utl_bilin_e
430
431
432	SUBROUTINE icb_utl_add( bergvals, ptvals )
433	!!----------------------------------------------------------------------
434	!! * ROUTINE icb_utl_add *
435	!!
436	!! ** Purpose : add a new berg to the iceberg list
437	!!
438	!!----------------------------------------------------------------------
439	TYPE(iceberg), INTENT(in) :: bergvals
440	TYPE(point) , INTENT(in) :: ptvals
441	!
442	TYPE(iceberg), POINTER :: new => NULL()
443	!!----------------------------------------------------------------------
444	!
445	new => NULL()
446	CALL icb_utl_create( new, bergvals, ptvals )
447	CALL icb_utl_insert( new )
448	new => NULL() ! Clear new
449	!
450	END SUBROUTINE icb_utl_add
451
452
453	SUBROUTINE icb_utl_create( berg, bergvals, ptvals )
454	!!----------------------------------------------------------------------
455	!! * ROUTINE icb_utl_create *
456	!!
457	!! ** Purpose : add a new berg to the iceberg list
458	!!
459	!!----------------------------------------------------------------------
460	TYPE(iceberg), INTENT(in) :: bergvals
461	TYPE(point) , INTENT(in) :: ptvals
462	TYPE(iceberg), POINTER :: berg
463	!
464	TYPE(point) , POINTER :: pt
465	INTEGER :: istat
466	!!----------------------------------------------------------------------
467	!
468	IF( ASSOCIATED(berg) ) CALL ctl_stop( 'icebergs, icb_utl_create: berg already associated' )
469	ALLOCATE(berg, STAT=istat)
470	IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate iceberg' )
471	berg%number(:) = bergvals%number(:)
472	berg%mass_scaling = bergvals%mass_scaling
473	berg%prev => NULL()
474	berg%next => NULL()
475	!
476	ALLOCATE(pt, STAT=istat)
477	IF( istat /= 0 ) CALL ctl_stop( 'failed to allocate first iceberg point' )
478	pt = ptvals
479	berg%current_point => pt
480	!
481	END SUBROUTINE icb_utl_create
482
483
484	SUBROUTINE icb_utl_insert( newberg )
485	!!----------------------------------------------------------------------
486	!! * ROUTINE icb_utl_insert *
487	!!
488	!! ** Purpose : add a new berg to the iceberg list
489	!!
490	!!----------------------------------------------------------------------
491	TYPE(iceberg), POINTER :: newberg
492	!
493	TYPE(iceberg), POINTER :: this, prev, last
494	!!----------------------------------------------------------------------
495	!
496	IF( ASSOCIATED( first_berg ) ) THEN
497	last => first_berg
498	DO WHILE (ASSOCIATED(last%next))
499	last => last%next
500	ENDDO
501	newberg%prev => last
502	last%next => newberg
503	last => newberg
504	ELSE ! list is empty so create it
505	first_berg => newberg
506	ENDIF
507	!
508	END SUBROUTINE icb_utl_insert
509
510
511	REAL(wp) FUNCTION icb_utl_yearday(kmon, kday, khr, kmin, ksec)
512	!!----------------------------------------------------------------------
513	!! * FUNCTION icb_utl_yearday *
514	!!
515	!! ** Purpose :
516	!!
517	! sga - improved but still only applies to 365 day year, need to do this properly
518	!
519	!!gm all these info are already known in daymod, no???
520	!!
521	!!----------------------------------------------------------------------
522	INTEGER, INTENT(in) :: kmon, kday, khr, kmin, ksec
523	!
524	INTEGER, DIMENSION(12) :: imonths = (/ 0,31,28,31,30,31,30,31,31,30,31,30 /)
525	!!----------------------------------------------------------------------
526	!
527	icb_utl_yearday = REAL( SUM( imonths(1:kmon) ), wp )
528	icb_utl_yearday = icb_utl_yearday + REAL(kday-1,wp) + (REAL(khr,wp) + (REAL(kmin,wp) + REAL(ksec,wp)/60.)/60.)/24.
529	!
530	END FUNCTION icb_utl_yearday
531
532	!!-------------------------------------------------------------------------
533
534	SUBROUTINE icb_utl_delete( first, berg )
535	!!----------------------------------------------------------------------
536	!! * ROUTINE icb_utl_delete *
537	!!
538	!! ** Purpose :
539	!!
540	!!----------------------------------------------------------------------
541	TYPE(iceberg), POINTER :: first, berg
542	!!----------------------------------------------------------------------
543	! Connect neighbors to each other
544	IF ( ASSOCIATED(berg%prev) ) THEN
545	berg%prev%next => berg%next
546	ELSE
547	first => berg%next
548	ENDIF
549	IF (ASSOCIATED(berg%next)) berg%next%prev => berg%prev
550	!
551	CALL icb_utl_destroy(berg)
552	!
553	END SUBROUTINE icb_utl_delete
554
555
556	SUBROUTINE icb_utl_destroy( berg )
557	!!----------------------------------------------------------------------
558	!! * ROUTINE icb_utl_destroy *
559	!!
560	!! ** Purpose : remove a single iceberg instance
561	!!
562	!!----------------------------------------------------------------------
563	TYPE(iceberg), POINTER :: berg
564	!!----------------------------------------------------------------------
565	!
566	! Remove any points
567	IF( ASSOCIATED( berg%current_point ) ) DEALLOCATE( berg%current_point )
568	!
569	DEALLOCATE(berg)
570	!
571	END SUBROUTINE icb_utl_destroy
572
573
574	SUBROUTINE icb_utl_track( knum, cd_label, kt )
575	!!----------------------------------------------------------------------
576	!! * ROUTINE icb_utl_track *
577	!!
578	!! ** Purpose :
579	!!
580	!!----------------------------------------------------------------------
581	INTEGER, DIMENSION(nkounts) :: knum ! iceberg number
582	CHARACTER(len=*) :: cd_label !
583	INTEGER :: kt ! timestep number
584	!
585	TYPE(iceberg), POINTER :: this
586	LOGICAL :: match
587	INTEGER :: k
588	!!----------------------------------------------------------------------
589	!
590	this => first_berg
591	DO WHILE( ASSOCIATED(this) )
592	match = .TRUE.
593	DO k = 1, nkounts
594	IF( this%number(k) /= knum(k) ) match = .FALSE.
595	END DO
596	IF( match ) CALL icb_utl_print_berg(this, kt)
597	this => this%next
598	END DO
599	!
600	END SUBROUTINE icb_utl_track
601
602
603	SUBROUTINE icb_utl_print_berg( berg, kt )
604	!!----------------------------------------------------------------------
605	!! * ROUTINE icb_utl_print_berg *
606	!!
607	!! ** Purpose : print one
608	!!
609	!!----------------------------------------------------------------------
610	TYPE(iceberg), POINTER :: berg
611	TYPE(point) , POINTER :: pt
612	INTEGER :: kt ! timestep number
613	!!----------------------------------------------------------------------
614	!
615	IF (nn_verbose_level == 0) RETURN
616	pt => berg%current_point
617	WRITE(numicb, 9200) kt, berg%number(1), &
618	pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel, &
619	pt%uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi
620	CALL flush( numicb )
621	9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4))
622	!
623	END SUBROUTINE icb_utl_print_berg
624
625
626	SUBROUTINE icb_utl_print( cd_label, kt )
627	!!----------------------------------------------------------------------
628	!! * ROUTINE icb_utl_print *
629	!!
630	!! ** Purpose : print many
631	!!
632	!!----------------------------------------------------------------------
633	CHARACTER(len=*) :: cd_label
634	INTEGER :: kt ! timestep number
635	!
636	INTEGER :: ibergs, inbergs
637	TYPE(iceberg), POINTER :: this
638	!!----------------------------------------------------------------------
639	!
640	IF (nn_verbose_level == 0) RETURN
641	this => first_berg
642	IF( ASSOCIATED(this) ) THEN
643	WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea
644	WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' ) &
645	& 'timestep', 'number', 'xi,yj','lon,lat','u,v','uo,vo','ua,va','ui,vi'
646	ENDIF
647	DO WHILE( ASSOCIATED(this) )
648	CALL icb_utl_print_berg(this, kt)
649	this => this%next
650	END DO
651	ibergs = icb_utl_count()
652	inbergs = ibergs
653	CALL mpp_sum('icbutl', inbergs)
654	IF( ibergs > 0 ) WRITE(numicb,'(a," there are",i5," bergs out of",i6," on PE ",i4)') &
655	& cd_label, ibergs, inbergs, narea
656	!
657	END SUBROUTINE icb_utl_print
658
659
660	SUBROUTINE icb_utl_incr()
661	!!----------------------------------------------------------------------
662	!! * ROUTINE icb_utl_incr *
663	!!
664	!! ** Purpose :
665	!!
666	! Small routine for coping with very large integer values labelling icebergs
667	! num_bergs is a array of integers
668	! the first member is incremented in steps of jpnij starting from narea
669	! this means each iceberg is labelled with a unique number
670	! when this gets to the maximum allowed integer the second and subsequent members are
671	! used to count how many times the member before cycles
672	!!----------------------------------------------------------------------
673	INTEGER :: ii, ibig
674	!!----------------------------------------------------------------------
675
676	ibig = HUGE(num_bergs(1))
677	IF( ibig-jpnij < num_bergs(1) ) THEN
678	num_bergs(1) = narea
679	DO ii = 2,nkounts
680	IF( num_bergs(ii) == ibig ) THEN
681	num_bergs(ii) = 0
682	IF( ii == nkounts ) CALL ctl_stop('Sorry, run out of iceberg number space')
683	ELSE
684	num_bergs(ii) = num_bergs(ii) + 1
685	EXIT
686	ENDIF
687	END DO
688	ELSE
689	num_bergs(1) = num_bergs(1) + jpnij
690	ENDIF
691	!
692	END SUBROUTINE icb_utl_incr
693
694
695	INTEGER FUNCTION icb_utl_count()
696	!!----------------------------------------------------------------------
697	!! * FUNCTION icb_utl_count *
698	!!
699	!! ** Purpose :
700	!!----------------------------------------------------------------------
701	TYPE(iceberg), POINTER :: this
702	!!----------------------------------------------------------------------
703	!
704	icb_utl_count = 0
705	this => first_berg
706	DO WHILE( ASSOCIATED(this) )
707	icb_utl_count = icb_utl_count+1
708	this => this%next
709	END DO
710	!
711	END FUNCTION icb_utl_count
712
713
714	REAL(wp) FUNCTION icb_utl_mass( first, justbits, justbergs )
715	!!----------------------------------------------------------------------
716	!! * FUNCTION icb_utl_mass *
717	!!
718	!! ** Purpose : compute the mass all iceberg, all berg bits or all bergs.
719	!!----------------------------------------------------------------------
720	TYPE(iceberg) , POINTER :: first
721	TYPE(point) , POINTER :: pt
722	LOGICAL, INTENT(in), OPTIONAL :: justbits, justbergs
723	!
724	TYPE(iceberg), POINTER :: this
725	!!----------------------------------------------------------------------
726	icb_utl_mass = 0._wp
727	this => first
728	!
729	IF( PRESENT( justbergs ) ) THEN
730	DO WHILE( ASSOCIATED( this ) )
731	pt => this%current_point
732	icb_utl_mass = icb_utl_mass + pt%mass * this%mass_scaling
733	this => this%next
734	END DO
735	ELSEIF( PRESENT(justbits) ) THEN
736	DO WHILE( ASSOCIATED( this ) )
737	pt => this%current_point
738	icb_utl_mass = icb_utl_mass + pt%mass_of_bits * this%mass_scaling
739	this => this%next
740	END DO
741	ELSE
742	DO WHILE( ASSOCIATED( this ) )
743	pt => this%current_point
744	icb_utl_mass = icb_utl_mass + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling
745	this => this%next
746	END DO
747	ENDIF
748	!
749	END FUNCTION icb_utl_mass
750
751
752	REAL(wp) FUNCTION icb_utl_heat( first, justbits, justbergs )
753	!!----------------------------------------------------------------------
754	!! * FUNCTION icb_utl_heat *
755	!!
756	!! ** Purpose : compute the heat in all iceberg, all bergies or all bergs.
757	!!----------------------------------------------------------------------
758	TYPE(iceberg) , POINTER :: first
759	LOGICAL, INTENT(in), OPTIONAL :: justbits, justbergs
760	!
761	TYPE(iceberg) , POINTER :: this
762	TYPE(point) , POINTER :: pt
763	!!----------------------------------------------------------------------
764	icb_utl_heat = 0._wp
765	this => first
766	!
767	IF( PRESENT( justbergs ) ) THEN
768	DO WHILE( ASSOCIATED( this ) )
769	pt => this%current_point
770	icb_utl_heat = icb_utl_heat + pt%mass * this%mass_scaling * pt%heat_density
771	this => this%next
772	END DO
773	ELSEIF( PRESENT(justbits) ) THEN
774	DO WHILE( ASSOCIATED( this ) )
775	pt => this%current_point
776	icb_utl_heat = icb_utl_heat + pt%mass_of_bits * this%mass_scaling * pt%heat_density
777	this => this%next
778	END DO
779	ELSE
780	DO WHILE( ASSOCIATED( this ) )
781	pt => this%current_point
782	icb_utl_heat = icb_utl_heat + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling * pt%heat_density
783	this => this%next
784	END DO
785	ENDIF
786	!
787	END FUNCTION icb_utl_heat
788
789	!!======================================================================
790	END MODULE icbutl

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

Download in other formats: