source: NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90 @ 13357

Last change on this file since 13357 was 13357, checked in by mathiot, 8 weeks ago

ticket #1900: add changes for Merino 2016 works. Results unchanged if flag ln_M2016 is set to .false. . Remaining step: test ln_M2016 = .true.

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