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

Last change on this file since 11738 was 11738, checked in by marc, 13 months ago

The Dr Hook changes from my perl code.

File size: 12.2 KB
Line 
1MODULE icbclv
2
3   !!======================================================================
4   !!                       ***  MODULE  icbclv  ***
5   !! Icebergs:  calving routines for iceberg calving
6   !!======================================================================
7   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
8   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
9   !!            -    !                            Removal of mapping from another grid
10   !!            -    !  2011-04  (Alderson)       Split into separate modules
11   !!            -    !  2011-05  (Alderson)       budgets into separate module
12   !!----------------------------------------------------------------------
13   !!----------------------------------------------------------------------
14   !!   icb_clv_flx   : transfer input flux of ice into iceberg classes
15   !!   icb_clv       : calve icebergs from stored ice
16   !!----------------------------------------------------------------------
17   USE par_oce        ! NEMO parameters
18   USE dom_oce        ! NEMO ocean domain
19   USE phycst         ! NEMO physical constants
20   USE lib_mpp        ! NEMO MPI library, lk_mpp in particular
21   USE lbclnk         ! NEMO boundary exchanges for gridded data
22
23   USE icb_oce        ! iceberg variables
24   USE icbdia         ! iceberg diagnostics
25   USE icbutl         ! iceberg utility routines
26
27   USE sbc_oce        ! for icesheet freshwater input variables
28   USE in_out_manager
29   USE iom
30
31   USE yomhook, ONLY: lhook, dr_hook
32   USE parkind1, ONLY: jprb, jpim
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   icb_clv_flx  ! routine called in icbstp.F90 module
38   PUBLIC   icb_clv      ! routine called in icbstp.F90 module
39
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE icb_clv_flx( kt )
48      !!----------------------------------------------------------------------
49      !!                 ***  ROUTINE icb_clv_flx  ***
50      !!
51      !! ** Purpose :   accumulate ice available for calving into class arrays
52      !!
53      !!----------------------------------------------------------------------
54      INTEGER, INTENT(in)             :: kt
55      !
56      REAL(wp)                        :: zcalving_used, zdist, zfact
57      REAL(wp)                        :: zgreenland_calving_sum, zantarctica_calving_sum
58      INTEGER                         :: jn, ji, jj                    ! loop counters
59      INTEGER                         :: imx                           ! temporary integer for max berg class
60      LOGICAL, SAVE                   :: ll_first_call = .TRUE.
61      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
62      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
63      REAL(KIND=jprb)               :: zhook_handle
64
65      CHARACTER(LEN=*), PARAMETER :: RoutineName='ICB_CLV_FLX'
66
67      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
68
69      !!----------------------------------------------------------------------
70      !
71      ! Adapt calving flux and calving heat flux from coupler for use here
72      ! Use interior mask: so no bergs in overlap areas and convert from km^3/year to kg/s
73      ! this assumes that input is given as equivalent water flux so that pure water density is appropriate
74
75      zfact = ( (1000._wp)**3 / ( NINT(rday) * nyear_len(1) ) ) * 850._wp
76      berg_grid%calving(:,:) = src_calving(:,:) * tmask_i(:,:) * zfact
77
78      IF( lk_oasis) THEN
79      ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
80      IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
81
82        ! Adjust total calving rates so that sum of iceberg calving and iceshelf melting in the northern
83        ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
84        ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
85
86         zgreenland_calving_sum = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) )
87         IF( lk_mpp ) CALL mpp_sum( zgreenland_calving_sum )
88         WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                                 &
89        &    berg_grid%calving(:,:) = berg_grid%calving(:,:) * greenland_icesheet_mass_rate_of_change * rn_greenland_calving_fraction &
90        &                                     / ( zgreenland_calving_sum + 1.0e-10_wp )
91
92         ! check
93         IF(lwp) WRITE(numout, *) 'Greenland iceberg calving climatology (kg/s) : ',zgreenland_calving_sum
94         zgreenland_calving_sum = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) )
95         IF( lk_mpp ) CALL mpp_sum( zgreenland_calving_sum )
96         IF(lwp) WRITE(numout, *) 'Greenland iceberg calving adjusted value (kg/s) : ',zgreenland_calving_sum
97
98         zantarctica_calving_sum = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) )
99         IF( lk_mpp ) CALL mpp_sum( zantarctica_calving_sum )
100         WHERE( antarctica_icesheet_mask(:,:) == 1.0 )                                                                              &
101         berg_grid%calving(:,:) = berg_grid%calving(:,:) * antarctica_icesheet_mass_rate_of_change * rn_antarctica_calving_fraction &
102        &                           / ( zantarctica_calving_sum + 1.0e-10_wp )
103 
104         ! check
105         IF(lwp) WRITE(numout, *) 'Antarctica iceberg calving climatology (kg/s) : ',zantarctica_calving_sum
106         zantarctica_calving_sum = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) )
107         IF( lk_mpp ) CALL mpp_sum( zantarctica_calving_sum )
108         IF(lwp) WRITE(numout, *) 'Antarctica iceberg calving adjusted value (kg/s) : ',zantarctica_calving_sum
109
110      ENDIF
111      ENDIF
112   
113      CALL iom_put( 'berg_calve', berg_grid%calving(:,:) )
114
115      ! Heat in units of W/m2, and mask (just in case)
116      berg_grid%calving_hflx(:,:) = src_calving_hflx(:,:) * tmask_i(:,:)
117
118      IF( ll_first_call .AND. .NOT. l_restarted_bergs) THEN      ! This is a hack to simplify initialization
119         ll_first_call = .FALSE.
120         !do jn=1, nclasses
121         !  where (berg_grid%calving==0.) berg_grid%stored_ice(:,:,jn)=0.
122         !end do
123         DO jj = 2, jpjm1
124            DO ji = 2, jpim1
125               IF( berg_grid%calving(ji,jj) /= 0._wp )                                  &    ! Need units of J
126                  berg_grid%stored_heat(ji,jj) = SUM( berg_grid%stored_ice(ji,jj,:) ) *         &  ! initial stored ice in kg
127                                         berg_grid%calving_hflx(ji,jj) * e1e2t(ji,jj) /   &  ! J/s/m2 x m^2 = J/s
128                                         berg_grid%calving(ji,jj)                            ! /calving in kg/s
129            END DO
130         END DO
131      ENDIF
132
133      ! assume that all calving flux must be distributed even if distribution array does not sum
134      ! to one - this may not be what is intended, but it's what you've got
135      DO jj = 1,jpj
136         DO ji = 1,jpi
137            imx = berg_grid%maxclass(ji,jj)
138            zdist = SUM( rn_distribution(1:nclasses) ) / SUM( rn_distribution(1:imx) )
139            DO jn = 1, imx
140               berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn) + &
141                                          berg_dt * berg_grid%calving(ji,jj) * rn_distribution(jn) * zdist
142            END DO
143         END DO
144      END DO
145
146      ! before changing the calving, save the amount we're about to use and do budget
147      zcalving_used = SUM( berg_grid%calving(:,:) )
148      berg_grid%tmp(:,:) = berg_dt * berg_grid%calving_hflx(:,:) * e1e2t(:,:) * tmask_i(:,:)
149      berg_grid%stored_heat (:,:) = berg_grid%stored_heat (:,:) + berg_grid%tmp(:,:)
150      CALL icb_dia_income( kt,  zcalving_used, berg_grid%tmp )
151      !
152      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
153   END SUBROUTINE icb_clv_flx
154
155   SUBROUTINE icb_clv()
156      !!----------------------------------------------------------------------
157      !!                 ***  ROUTINE icb_clv  ***
158      !!
159      !! ** Purpose :   This routine takes a stored ice field and calves to the ocean,
160      !!                so the gridded array stored_ice has only non-zero entries at selected
161      !!                wet points adjacent to known land based calving points
162      !!
163      !! ** method  : - Look at each grid point and see if there's enough for each size class to calve
164      !!                If there is, a new iceberg is calved.  This happens in the order determined by
165      !!                the class definition arrays (which in the default case is smallest first)
166      !!                Note that only the non-overlapping part of the processor where icebergs are allowed
167      !!                is considered
168      !!----------------------------------------------------------------------
169      INTEGER       ::   ji, jj, jn   ! dummy loop indices
170      INTEGER       ::   icnt, icntmax
171      TYPE(iceberg) ::   newberg
172      TYPE(point)   ::   newpt
173      REAL(wp)      ::   zday, zcalved_to_berg, zheat_to_berg
174      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
175      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
176      REAL(KIND=jprb)               :: zhook_handle
177
178      CHARACTER(LEN=*), PARAMETER :: RoutineName='ICB_CLV'
179
180      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
181
182      !!----------------------------------------------------------------------
183      !
184      icntmax = 0
185      zday    = REAL(nday_year,wp) + REAL(nsec_day,wp)/86400.0_wp
186      !
187      DO jn = 1, nclasses
188         DO jj = nicbdj, nicbej
189            DO ji = nicbdi, nicbei
190               !
191               icnt = 0
192               !
193               DO WHILE (berg_grid%stored_ice(ji,jj,jn) >= rn_initial_mass(jn) * rn_mass_scaling(jn) )
194                  !
195                  newpt%lon = glamt(ji,jj)         ! at t-point (centre of the cell)
196                  newpt%lat = gphit(ji,jj)
197                  newpt%xi  = REAL( mig(ji), wp )
198                  newpt%yj  = REAL( mjg(jj), wp )
199                  !
200                  newpt%uvel = 0._wp               ! initially at rest
201                  newpt%vvel = 0._wp
202                  !                                ! set berg characteristics
203                  newpt%mass           = rn_initial_mass     (jn)
204                  newpt%thickness      = rn_initial_thickness(jn)
205                  newpt%width          = first_width         (jn)
206                  newpt%length         = first_length        (jn)
207                  newberg%mass_scaling = rn_mass_scaling     (jn)
208                  newpt%mass_of_bits   = 0._wp                          ! no bergy
209                  !
210                  newpt%year   = nyear
211                  newpt%day    = zday
212                  newpt%heat_density = berg_grid%stored_heat(ji,jj) / berg_grid%stored_ice(ji,jj,jn)   ! This is in J/kg
213                  !
214                  CALL icb_utl_incr()
215                  newberg%number(:) = num_bergs(:)
216                  !
217                  CALL icb_utl_add( newberg, newpt )
218                  !
219                  zcalved_to_berg = rn_initial_mass(jn) * rn_mass_scaling(jn)           ! Units of kg
220                  !                                ! Heat content
221                  zheat_to_berg           = zcalved_to_berg * newpt%heat_density             ! Units of J
222                  berg_grid%stored_heat(ji,jj) = berg_grid%stored_heat(ji,jj) - zheat_to_berg
223                  !                                ! Stored mass
224                  berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn) - zcalved_to_berg
225                  !
226                  icnt = icnt + 1
227                  !
228                  CALL icb_dia_calve(ji, jj, jn,  zcalved_to_berg, zheat_to_berg )
229               END DO
230               icntmax = MAX( icntmax, icnt )
231            END DO
232         END DO
233      END DO
234      !
235      DO jn = 1,nclasses
236         CALL lbc_lnk( berg_grid%stored_ice(:,:,jn), 'T', 1._wp )
237      END DO
238      CALL lbc_lnk( berg_grid%stored_heat, 'T', 1._wp )
239      !
240      IF( nn_verbose_level > 0 .AND. icntmax > 1 )   WRITE(numicb,*) 'icb_clv: icnt=', icnt,' on', narea
241      !
242      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
243   END SUBROUTINE  icb_clv
244
245   !!======================================================================
246END MODULE icbclv
Note: See TracBrowser for help on using the repository browser.