source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/ICB/icbutl.F90 @ 10402

Last change on this file since 10402 was 10402, checked in by smasson, 21 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: no more need of lk_mpp for mpp_sum/max/min, see #2133

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