Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

icbutl.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago
The Dr Hook changes from my perl code.
File size: 38.4 KB

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90 @ 11738

Download in other formats: