source: NEMO/trunk/src/OCE/ICB/icbutl.F90 @ 10570

Last change on this file since 10570 was 10570, checked in by acc, 19 months ago

Trunk update to implement finer control over the choice of text report files generated. See ticket: #2167

  • 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      IF (nn_verbose_level == 0) RETURN
628      pt => berg%current_point
629      WRITE(numicb, 9200) kt, berg%number(1), &
630                   pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel,  &
631                   pt%uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi
632      CALL flush( numicb )
633 9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4))
634      !
635   END SUBROUTINE icb_utl_print_berg
636
637
638   SUBROUTINE icb_utl_print( cd_label, kt )
639      !!----------------------------------------------------------------------
640      !!                 ***  ROUTINE icb_utl_print  ***
641      !!
642      !! ** Purpose :   print many
643      !!
644      !!----------------------------------------------------------------------
645      CHARACTER(len=*)       :: cd_label
646      INTEGER                :: kt             ! timestep number
647      !
648      INTEGER                :: ibergs, inbergs
649      TYPE(iceberg), POINTER :: this
650      !!----------------------------------------------------------------------
651      !
652      IF (nn_verbose_level == 0) RETURN
653      this => first_berg
654      IF( ASSOCIATED(this) ) THEN
655         WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea
656         WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' )   &
657            &         'timestep', 'number', 'xi,yj','lon,lat','u,v','uo,vo','ua,va','ui,vi'
658      ENDIF
659      DO WHILE( ASSOCIATED(this) )
660        CALL icb_utl_print_berg(this, kt)
661        this => this%next
662      END DO
663      ibergs = icb_utl_count()
664      inbergs = ibergs
665      CALL mpp_sum('icbutl', inbergs)
666      IF( ibergs > 0 )   WRITE(numicb,'(a," there are",i5," bergs out of",i6," on PE ",i4)')   &
667         &                                  cd_label, ibergs, inbergs, narea
668      !
669   END SUBROUTINE icb_utl_print
670
671
672   SUBROUTINE icb_utl_incr()
673      !!----------------------------------------------------------------------
674      !!                 ***  ROUTINE icb_utl_incr  ***
675      !!
676      !! ** Purpose :   
677      !!
678      ! Small routine for coping with very large integer values labelling icebergs
679      ! num_bergs is a array of integers
680      ! the first member is incremented in steps of jpnij starting from narea
681      ! this means each iceberg is labelled with a unique number
682      ! when this gets to the maximum allowed integer the second and subsequent members are
683      ! used to count how many times the member before cycles
684      !!----------------------------------------------------------------------
685      INTEGER ::   ii, ibig
686      !!----------------------------------------------------------------------
687
688      ibig = HUGE(num_bergs(1))
689      IF( ibig-jpnij < num_bergs(1) ) THEN
690         num_bergs(1) = narea
691         DO ii = 2,nkounts
692            IF( num_bergs(ii) == ibig ) THEN
693               num_bergs(ii) = 0
694               IF( ii == nkounts ) CALL ctl_stop('Sorry, run out of iceberg number space')
695            ELSE
696               num_bergs(ii) = num_bergs(ii) + 1
697               EXIT
698            ENDIF
699         END DO
700      ELSE
701         num_bergs(1) = num_bergs(1) + jpnij
702      ENDIF
703      !
704   END SUBROUTINE icb_utl_incr
705
706
707   INTEGER FUNCTION icb_utl_count()
708      !!----------------------------------------------------------------------
709      !!                 ***  FUNCTION icb_utl_count  ***
710      !!
711      !! ** Purpose :   
712      !!----------------------------------------------------------------------
713      TYPE(iceberg), POINTER :: this
714      !!----------------------------------------------------------------------
715      !
716      icb_utl_count = 0
717      this => first_berg
718      DO WHILE( ASSOCIATED(this) )
719         icb_utl_count = icb_utl_count+1
720         this => this%next
721      END DO
722      !
723   END FUNCTION icb_utl_count
724
725
726   REAL(wp) FUNCTION icb_utl_mass( first, justbits, justbergs )
727      !!----------------------------------------------------------------------
728      !!                 ***  FUNCTION icb_utl_mass  ***
729      !!
730      !! ** Purpose :   compute the mass all iceberg, all berg bits or all bergs.
731      !!----------------------------------------------------------------------
732      TYPE(iceberg)      , POINTER  ::   first
733      TYPE(point)        , POINTER  ::   pt
734      LOGICAL, INTENT(in), OPTIONAL ::   justbits, justbergs
735      !
736      TYPE(iceberg), POINTER ::   this
737      !!----------------------------------------------------------------------
738      icb_utl_mass = 0._wp
739      this => first
740      !
741      IF( PRESENT( justbergs  ) ) THEN
742         DO WHILE( ASSOCIATED( this ) )
743            pt => this%current_point
744            icb_utl_mass = icb_utl_mass + pt%mass         * this%mass_scaling
745            this => this%next
746         END DO
747      ELSEIF( PRESENT(justbits) ) THEN
748         DO WHILE( ASSOCIATED( this ) )
749            pt => this%current_point
750            icb_utl_mass = icb_utl_mass + pt%mass_of_bits * this%mass_scaling
751            this => this%next
752         END DO
753      ELSE
754         DO WHILE( ASSOCIATED( this ) )
755            pt => this%current_point
756            icb_utl_mass = icb_utl_mass + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling
757            this => this%next
758         END DO
759      ENDIF
760      !
761   END FUNCTION icb_utl_mass
762
763
764   REAL(wp) FUNCTION icb_utl_heat( first, justbits, justbergs )
765      !!----------------------------------------------------------------------
766      !!                 ***  FUNCTION icb_utl_heat  ***
767      !!
768      !! ** Purpose :   compute the heat in all iceberg, all bergies or all bergs.
769      !!----------------------------------------------------------------------
770      TYPE(iceberg)      , POINTER  ::   first
771      LOGICAL, INTENT(in), OPTIONAL ::   justbits, justbergs
772      !
773      TYPE(iceberg)      , POINTER  ::   this
774      TYPE(point)        , POINTER  ::   pt
775      !!----------------------------------------------------------------------
776      icb_utl_heat = 0._wp
777      this => first
778      !
779      IF( PRESENT( justbergs  ) ) THEN
780         DO WHILE( ASSOCIATED( this ) )
781            pt => this%current_point
782            icb_utl_heat = icb_utl_heat + pt%mass         * this%mass_scaling * pt%heat_density
783            this => this%next
784         END DO
785      ELSEIF( PRESENT(justbits) ) THEN
786         DO WHILE( ASSOCIATED( this ) )
787            pt => this%current_point
788            icb_utl_heat = icb_utl_heat + pt%mass_of_bits * this%mass_scaling * pt%heat_density
789            this => this%next
790         END DO
791      ELSE
792         DO WHILE( ASSOCIATED( this ) )
793            pt => this%current_point
794            icb_utl_heat = icb_utl_heat + ( pt%mass + pt%mass_of_bits ) * this%mass_scaling * pt%heat_density
795            this => this%next
796         END DO
797      ENDIF
798      !
799   END FUNCTION icb_utl_heat
800
801   !!======================================================================
802END MODULE icbutl
Note: See TracBrowser for help on using the repository browser.