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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
iceitd.F90 in branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90 @ 8892

Last change on this file since 8892 was 8892, checked in by frrh, 7 years ago

Commit updates with debugging write statements.
File size: 38.5 KB
Line 
1MODULE iceitd
2   !!======================================================================
3   !!                       ***  MODULE iceitd ***
4   !!   sea-ice : ice thickness distribution
5   !!======================================================================
6   !! History :   -   !          (W. H. Lipscomb and E.C. Hunke) CICE (c) original code
7   !!            3.0  ! 2005-12  (M. Vancoppenolle) adaptation to LIM-3
8   !!             -   ! 2006-06  (M. Vancoppenolle) adaptation to include salt, age
9   !!             -   ! 2007-04  (M. Vancoppenolle) Mass conservation checked
10   !!----------------------------------------------------------------------
11#if defined key_lim3
12   !!----------------------------------------------------------------------
13   !!   'key_lim3'                                       ESIM sea-ice model
14   !!----------------------------------------------------------------------
15   !!   ice_itd_init   : read ice thicknesses mean and min from namelist
16   !!   ice_itd_rem    : redistribute ice thicknesses after thermo growth and melt
17   !!   ice_itd_reb    : rebin ice thicknesses into bounded categories
18   !!----------------------------------------------------------------------
19   USE dom_oce        ! ocean domain
20   USE phycst         ! physical constants
21   USE ice1D          ! sea-ice: thermodynamic variables
22   USE ice            ! sea-ice: variables
23   USE icectl         ! sea-ice: conservation tests
24   USE icetab         ! sea-ice: convert 1D<=>2D
25   !
26   USE in_out_manager ! I/O manager
27   USE lib_mpp        ! MPP library
28   USE lib_fortran    ! fortran utilities (sum + no signed zero)
29   USE prtctl         ! Print control
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   ice_itd_init  ! called in icestp
35   PUBLIC   ice_itd_rem   ! called in icethd
36   PUBLIC   ice_itd_reb   ! called in icecor
37
38   INTEGER ::              nice_catbnd   ! choice of the type of ice category function
39   !                                       ! associated indices:
40   INTEGER, PARAMETER ::   np_cathfn = 1   ! categories defined by a function
41   INTEGER, PARAMETER ::   np_catusr = 2   ! categories defined by the user
42   !
43   ! ** namelist (namitd) **
44   LOGICAL  ::   ln_cat_hfn  ! ice categories are defined by function like rn_himean**(-0.05)
45   REAL(wp) ::   rn_himean   ! mean thickness of the domain
46   LOGICAL  ::   ln_cat_usr  ! ice categories are defined by rn_catbnd
47   REAL(wp), DIMENSION(0:100) ::   rn_catbnd   ! ice categories bounds
48   !
49   !!----------------------------------------------------------------------
50   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
51   !! $Id: iceitd.F90 8420 2017-08-08 12:18:46Z clem $
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE ice_itd_rem( kt )
57      !!------------------------------------------------------------------
58      !!                ***  ROUTINE ice_itd_rem ***
59      !!
60      !! ** Purpose :   computes the redistribution of ice thickness
61      !!                after thermodynamic growth of ice thickness
62      !!
63      !! ** Method  :   Linear remapping
64      !!
65      !! References :   W.H. Lipscomb, JGR 2001
66      !!------------------------------------------------------------------
67      INTEGER , INTENT (in) ::   kt      ! Ocean time step
68      !
69      INTEGER  ::   ji, jj, jl, jcat     ! dummy loop index
70      INTEGER  ::   ipti                 ! local integer
71      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars
72      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
73      REAL(wp) ::   zx3       
74      REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds"
75
76      INTEGER , DIMENSION(jpij)       ::   iptidx          ! compute remapping or not
77      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor          ! donor category index
78      REAL(wp), DIMENSION(jpij,jpl)   ::   zdhice          ! ice thickness increment
79      REAL(wp), DIMENSION(jpij,jpl)   ::   g0, g1          ! coefficients for fitting the line of the ITD
80      REAL(wp), DIMENSION(jpij,jpl)   ::   hL, hR          ! left and right boundary for the ITD for each thickness
81      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice  ! local increment of ice area and volume
82      REAL(wp), DIMENSION(jpij)       ::   zhb0, zhb1      ! category boundaries for thinnes categories
83      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories
84      !!------------------------------------------------------------------
85
86      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' 
87
88      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
89
90      !-----------------------------------------------------------------------------------------------
91      !  1) Identify grid cells with ice
92      !-----------------------------------------------------------------------------------------------
93      npti = 0   ;   nptidx(:) = 0
94      DO jj = 1, jpj
95         DO ji = 1, jpi
96            IF ( at_i(ji,jj) > epsi10 ) THEN
97               npti         = npti + 1
98               nptidx( npti ) = (jj - 1) * jpi + ji
99            ENDIF
100         END DO
101      END DO
102     
103      !-----------------------------------------------------------------------------------------------
104      !  2) Compute new category boundaries
105      !-----------------------------------------------------------------------------------------------
106      IF( npti > 0 ) THEN
107
108         zdhice(:,:) = 0._wp
109         zhbnew(:,:) = 0._wp
110write(numout,*) "RSRH PE iceitd AA h_i", sum(h_i(:,:,:)); flush(numout)
111write(numout,*) "RSRH PE iceitd AA a_i", sum(a_i(:,:,:)); flush(numout)
112         CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i   )
113         CALL tab_3d_2d( npti, nptidx(1:npti), h_ib_2d(1:npti,1:jpl), h_i_b )
114         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i    )
115         CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d (1:npti,1:jpl), a_i_b  )
116
117         DO jl = 1, jpl
118            ! Compute thickness change in each ice category
119            DO ji = 1, npti
120               zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl)
121            END DO
122         END DO
123         
124         ! --- New boundaries for category 1:jpl-1 --- !
125         DO jl = 1, jpl - 1
126            !
127            DO ji = 1, npti
128               !
129               ! --- New boundary: Hn* = Hn + Fn*dt --- !
130               !     Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - h_i_b)
131               !
132               IF    ( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl+1) & a(jl) /= 0
133                  zslope           = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( h_ib_2d(ji,jl+1) - h_ib_2d(ji,jl) )
134                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - h_ib_2d(ji,jl) )
135               ELSEIF( a_ib_2d(ji,jl) >  epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN   ! a(jl+1)=0 => Hn* = Hn + fn*dt
136                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl)
137               ELSEIF( a_ib_2d(ji,jl) <= epsi10 .AND. a_ib_2d(ji,jl+1) >  epsi10 ) THEN   ! a(jl)=0 => Hn* = Hn + fn+1*dt
138                  zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl+1)
139               ELSE                                                                       ! a(jl+1) & a(jl) = 0
140                  zhbnew(ji,jl) = hi_max(jl)
141               ENDIF
142               !
143               ! --- 2 conditions for remapping --- !
144               ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi               
145               !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
146               !          in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
147               IF( a_i_2d(ji,jl  ) > epsi10 .AND. h_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi10 ) )   nptidx(ji) = 0
148               IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) )   nptidx(ji) = 0
149               
150               ! 2) Hn-1 < Hn* < Hn+1 
151               IF( zhbnew(ji,jl) < hi_max(jl-1) )   nptidx(ji) = 0
152               IF( zhbnew(ji,jl) > hi_max(jl+1) )   nptidx(ji) = 0
153               
154            END DO
155         END DO
156         !
157         ! --- New boundaries for category jpl --- !
158         DO ji = 1, npti
159            IF( a_i_2d(ji,jpl) > epsi10 ) THEN
160               zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * h_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) )
161            ELSE
162               zhbnew(ji,jpl) = hi_max(jpl) 
163            ENDIF
164           
165            ! --- 1 additional condition for remapping (1st category) --- !
166            ! H0+epsi < h1(t) < H1-epsi
167            !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
168            !    in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
169            IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) )   nptidx(ji) = 0
170            IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   nptidx(ji) = 0
171         END DO
172         !
173         !-----------------------------------------------------------------------------------------------
174         !  3) Identify cells where remapping
175         !-----------------------------------------------------------------------------------------------
176         ipti = 0 ; iptidx(:) = 0
177         DO ji = 1, npti
178            IF( nptidx(ji) /= 0 ) THEN
179               ipti = ipti + 1
180               iptidx(ipti)   = nptidx(ji)
181               zhbnew(ipti,:) = zhbnew(ji,:)  ! adjust zhbnew to new indices
182            ENDIF
183         END DO
184         nptidx(:) = iptidx(:)
185         npti      = ipti
186         !
187      ENDIF
188   
189      !-----------------------------------------------------------------------------------------------
190      !  4) Compute g(h)
191      !-----------------------------------------------------------------------------------------------
192      IF( npti > 0 ) THEN
193         !
194         zhb0(:) = hi_max(0)   ;   zhb1(:) = hi_max(1)
195         g0(:,:) = 0._wp       ;   g1(:,:) = 0._wp 
196         hL(:,:) = 0._wp       ;   hR(:,:) = 0._wp 
197         !
198         DO jl = 1, jpl
199            !
200write(numout,*) "RSRH PE iceitd BB h_i", sum(h_i(:,:,:)); flush(numout)
201write(numout,*) "RSRH PE iceitd BB a_i", sum(a_i(:,:,:)); flush(numout)
202            CALL tab_2d_1d( npti, nptidx(1:npti), h_ib_1d(1:npti), h_i_b(:,:,jl) )
203            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)   )
204            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    )
205            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl)    )
206            !
207            IF( jl == 1 ) THEN
208               ! 
209               ! --- g(h) for category 1 --- !
210               CALL ice_itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in
211                  &                  g0  (1:npti,1), g1  (1:npti,1), hL      (1:npti,1), hR    (1:npti,1)   )   ! out
212                  !
213               ! Area lost due to melting of thin ice
214               DO ji = 1, npti
215                  !
216                  IF( a_i_1d(ji) > epsi10 ) THEN
217                     !
218                     zdh0 =  h_i_1d(ji) - h_ib_1d(ji)               
219                     IF( zdh0 < 0.0 ) THEN      !remove area from category 1
220                        zdh0 = MIN( -zdh0, hi_max(1) )
221                        !Integrate g(1) from 0 to dh0 to estimate area melted
222                        zetamax = MIN( zdh0, hR(ji,1) ) - hL(ji,1)
223                        !
224                        IF( zetamax > 0.0 ) THEN
225                           zx1    = zetamax
226                           zx2    = 0.5 * zetamax * zetamax 
227                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                        ! ice area removed
228                           zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i               
229                           zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting
230                           !     of thin ice (zdamax > 0)
231                           ! Remove area, conserving volume
232                           h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )
233                           a_i_1d(ji)  = a_i_1d(ji) - zda0
234                           v_i_1d(ji)  = a_i_1d(ji) * h_i_1d(ji) ! clem-useless ?
235                        ENDIF
236                        !
237                     ELSE ! if ice accretion zdh0 > 0
238                        ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
239                        zhbnew(ji,0) = MIN( zdh0, hi_max(1) ) 
240                     ENDIF
241                     !
242                  ENDIF
243                  !
244               END DO
245               !
246write(numout,*) "RSRH PE iceitd CC h_i", sum(h_i(:,:,:)); flush(numout)
247write(numout,*) "RSRH PE iceitd CC a_i", sum(a_i(:,:,:)); flush(numout)
248               CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)   )
249               CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    )
250               CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl)    )
251               !
252            ENDIF ! jl=1
253            !
254            ! --- g(h) for each thickness category --- ! 
255            CALL ice_itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in
256               &                  g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR    (1:npti,jl)   )   ! out
257            !
258         END DO
259         
260         !-----------------------------------------------------------------------------------------------
261         !  5) Compute area and volume to be shifted across each boundary (Eq. 18)
262         !-----------------------------------------------------------------------------------------------
263         DO jl = 1, jpl - 1
264            !
265            DO ji = 1, npti
266               !
267               ! left and right integration limits in eta space
268               IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1
269                  zetamin = MAX( hi_max(jl)   , hL(ji,jl) ) - hL(ji,jl)   ! hi_max(jl) - hL
270                  zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl)   ! hR - hL
271                  jdonor(ji,jl) = jl
272               ELSE                                 ! Hn* <= Hn => transfer from jl+1 to jl
273                  zetamin = 0.0
274                  zetamax = MIN( hi_max(jl), hR(ji,jl+1) ) - hL(ji,jl+1)  ! hi_max(jl) - hL
275                  jdonor(ji,jl) = jl + 1
276               ENDIF
277               zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin
278               !
279               zx1  = zetamax - zetamin
280               zwk1 = zetamin * zetamin
281               zwk2 = zetamax * zetamax
282               zx2  = 0.5 * ( zwk2 - zwk1 )
283               zwk1 = zwk1 * zetamin
284               zwk2 = zwk2 * zetamax
285               zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 )
286               jcat = jdonor(ji,jl)
287               zdaice(ji,jl) = g1(ji,jcat)*zx2 + g0(ji,jcat)*zx1
288               zdvice(ji,jl) = g1(ji,jcat)*zx3 + g0(ji,jcat)*zx2 + zdaice(ji,jl)*hL(ji,jcat)
289               !
290            END DO
291         END DO
292         
293         !----------------------------------------------------------------------------------------------
294         ! 6) Shift ice between categories
295         !----------------------------------------------------------------------------------------------
296         CALL ice_itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )
297         
298         !----------------------------------------------------------------------------------------------
299         ! 7) Make sure h_i >= minimum ice thickness hi_min
300         !----------------------------------------------------------------------------------------------
301write(numout,*) "RSRH PE iceitd DD h_i", sum(h_i(:,:,:)); flush(numout)
302write(numout,*) "RSRH PE iceitd DD a_i", sum(a_i(:,:,:)); flush(numout)
303write(numout,*) "RSRH PE iceitd DD a_ip", sum(a_ip(:,:,:)); flush(numout)
304         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,1)   )
305         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,1)    )
306         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1)   )
307write(numout,*) "RSRH PE iceitd EE h_i", sum(h_i(:,:,:)); flush(numout)
308write(numout,*) "RSRH PE iceitd EE a_i", sum(a_i(:,:,:)); flush(numout)
309write(numout,*) "RSRH PE iceitd EE a_ip", sum(a_ip(:,:,:)); flush(numout)
310         
311         DO ji = 1, npti
312            IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
313               a_i_1d (ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 
314               IF ( ln_pnd_H12 )    a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
315               h_i_1d(ji) = rn_himin
316            ENDIF
317         END DO
318         !
319write(numout,*) "RSRH PE iceitd FF h_i", sum(h_i(:,:,:)); flush(numout)
320write(numout,*) "RSRH PE iceitd FF a_i", sum(a_i(:,:,:)); flush(numout)
321write(numout,*) "RSRH PE iceitd FF a_ip", sum(a_ip(:,:,:)); flush(numout)
322         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,1)   )
323         CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,1)    )
324         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1)   )
325write(numout,*) "RSRH PE iceitd GG h_i", sum(h_i(:,:,:)); flush(numout)
326write(numout,*) "RSRH PE iceitd GG a_i", sum(a_i(:,:,:)); flush(numout)
327write(numout,*) "RSRH PE iceitd GG a_ip", sum(a_ip(:,:,:)); flush(numout)
328         !
329      ENDIF
330      !
331      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
332      !
333   END SUBROUTINE ice_itd_rem
334
335
336   SUBROUTINE ice_itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR )
337      !!------------------------------------------------------------------
338      !!                ***  ROUTINE ice_itd_glinear ***
339      !!
340      !! ** Purpose :   build g(h) satisfying area and volume constraints (Eq. 6 and 9)
341      !!
342      !! ** Method  :   g(h) is linear and written as: g(eta) = g1(eta) + g0
343      !!                with eta = h - HL
344      !!------------------------------------------------------------------
345      REAL(wp), DIMENSION(:), INTENT(in   ) ::   HbL, HbR      ! left and right category boundaries
346      REAL(wp), DIMENSION(:), INTENT(in   ) ::   phice, paice  ! ice thickness and concentration
347      REAL(wp), DIMENSION(:), INTENT(inout) ::   pg0, pg1      ! coefficients in linear equation for g(eta)
348      REAL(wp), DIMENSION(:), INTENT(inout) ::   phL, phR      ! min and max value of range over which g(h) > 0
349      !
350      INTEGER  ::   ji           ! horizontal indices
351      REAL(wp) ::   z1_3 , z2_3  ! 1/3 , 2/3
352      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL)
353      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL)
354      REAL(wp) ::   zdhr         ! 1 / (hR - hL)
355      REAL(wp) ::   zwk1, zwk2   ! temporary variables
356      !!------------------------------------------------------------------
357      !
358      z1_3 = 1._wp / 3._wp
359      z2_3 = 2._wp / 3._wp
360      !
361      DO ji = 1, npti
362         !
363         IF( paice(ji) > epsi10  .AND. phice(ji) > 0._wp )  THEN
364            !
365            ! Initialize hL and hR
366            phL(ji) = HbL(ji)
367            phR(ji) = HbR(ji)
368            !
369            ! Change hL or hR if hice falls outside central third of range,
370            ! so that hice is in the central third of the range [HL HR]
371            zh13 = z1_3 * ( 2._wp * phL(ji) +         phR(ji) )
372            zh23 = z1_3 * (         phL(ji) + 2._wp * phR(ji) )
373            !
374            IF    ( phice(ji) < zh13 ) THEN   ;   phR(ji) = 3._wp * phice(ji) - 2._wp * phL(ji) ! move HR to the left
375            ELSEIF( phice(ji) > zh23 ) THEN   ;   phL(ji) = 3._wp * phice(ji) - 2._wp * phR(ji) ! move HL to the right
376            ENDIF
377            !
378            ! Compute coefficients of g(eta) = g0 + g1*eta
379            zdhr = 1._wp / (phR(ji) - phL(ji))
380            zwk1 = 6._wp * paice(ji) * zdhr
381            zwk2 = ( phice(ji) - phL(ji) ) * zdhr
382            pg0(ji) = zwk1 * ( z2_3 - zwk2 )                    ! Eq. 14
383            pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5_wp )   ! Eq. 14
384            !
385         ELSE  ! remap_flag = .false. or a_i < epsi10
386            phL(ji) = 0._wp
387            phR(ji) = 0._wp
388            pg0(ji) = 0._wp
389            pg1(ji) = 0._wp
390         ENDIF
391         !
392      END DO
393      !
394   END SUBROUTINE ice_itd_glinear
395
396
397   SUBROUTINE ice_itd_shiftice( kdonor, pdaice, pdvice )
398      !!------------------------------------------------------------------
399      !!                ***  ROUTINE ice_itd_shiftice ***
400      !!
401      !! ** Purpose :   shift ice across category boundaries, conserving everything
402      !!              ( area, volume, energy, age*vol, and mass of salt )
403      !!------------------------------------------------------------------
404      INTEGER , DIMENSION(:,:), INTENT(in) ::   kdonor   ! donor category index
405      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdaice   ! ice area transferred across boundary
406      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary
407      !
408      INTEGER  ::   ji, jj, jl, jk     ! dummy loop indices
409      INTEGER  ::   ii, ij, jl2, jl1   ! local integers
410      REAL(wp) ::   ztrans             ! ice/snow transferred
411      REAL(wp), DIMENSION(jpij)     ::   zworka, zworkv   ! workspace
412      REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn           !  -    -
413      !!------------------------------------------------------------------
414         
415      CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i  )
416      CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i  )
417      CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i  )
418      CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s  )
419      CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i )
420      CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i )
421      CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
422      CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
423      CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
424
425      !----------------------------------------------------------------------------------------------
426      ! 1) Define a variable equal to a_i*T_su
427      !----------------------------------------------------------------------------------------------
428      DO jl = 1, jpl
429         DO ji = 1, npti
430            zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl)
431         END DO
432      END DO
433     
434      !-------------------------------------------------------------------------------
435      ! 2) Transfer volume and energy between categories
436      !-------------------------------------------------------------------------------
437      DO jl = 1, jpl - 1
438         DO ji = 1, npti
439            !
440            jl1 = kdonor(ji,jl)
441            !
442            IF( jl1 > 0 ) THEN
443               !
444               IF ( jl1 == jl  ) THEN   ;   jl2 = jl1+1
445               ELSE                     ;   jl2 = jl 
446               ENDIF
447               !
448               IF( v_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
449               ELSE                                  ;   zworkv(ji) = 0._wp
450               ENDIF
451               IF( a_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
452               ELSE                                  ;   zworka(ji) = 0._wp
453               ENDIF
454               !
455               a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
456               a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
457               !
458               v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)       ! Ice volumes
459               v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
460               !
461               ztrans         = v_s_2d(ji,jl1) * zworkv(ji)          ! Snow volumes
462               v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
463               v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 
464               !         
465               !                                                     ! Ice age
466               ztrans          = oa_i_2d(ji,jl1) * pdaice(ji,jl)  !!clem: should be * zworka(ji) but it does not work ????
467               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - ztrans
468               oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ztrans
469               !
470               ztrans          = sv_i_2d(ji,jl1) * zworkv(ji)        ! Ice salinity
471               !
472               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - ztrans
473               sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ztrans
474               !
475               !                                                     ! Surface temperature
476               ztrans          = t_su_2d(ji,jl1) * pdaice(ji,jl)     !!clem: should be * zworka(ji) but it does not work ????
477               zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
478               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
479               ! 
480               IF ( ln_pnd_H12 ) THEN
481                  !                                                  ! Pond fraction
482                  ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
483                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
484                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
485                  !                                                  ! Pond volume (also proportional to da/a)
486                  ztrans          = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
487                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
488                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
489               ENDIF
490               !
491            ENDIF   ! jl1 >0
492         END DO
493         !
494         DO jk = 1, nlay_s         !--- Snow heat content
495            !
496            DO ji = 1, npti
497               ii = MOD( nptidx(ji) - 1, jpi ) + 1
498               ij = ( nptidx(ji) - 1 ) / jpi + 1
499               !
500               jl1 = kdonor(ji,jl)
501               !
502               IF( jl1 > 0 ) THEN
503                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
504                  ELSE                ;  jl2 = jl
505                  ENDIF
506                  !
507                  ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
508                  e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
509                  e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
510               ENDIF
511            END DO
512         END DO
513
514         DO jk = 1, nlay_i         !--- Ice heat content
515            DO ji = 1, npti
516               ii = MOD( nptidx(ji) - 1, jpi ) + 1
517               ij = ( nptidx(ji) - 1 ) / jpi + 1
518               !
519               jl1 = kdonor(ji,jl)
520               !
521               IF( jl1 > 0 ) THEN
522                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
523                  ELSE                ;  jl2 = jl
524                  ENDIF
525                  !
526                  ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji)
527                  e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
528                  e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
529               ENDIF
530            END DO
531         END DO
532         !
533      END DO                   ! boundaries, 1 to jpl-1
534     
535      !-------------------------------------------------------------------------------
536      ! 3) Update ice thickness and temperature
537      !-------------------------------------------------------------------------------
538      WHERE( a_i_2d(1:npti,:) >= epsi20 )
539         h_i_2d(1:npti,:)  =  v_i_2d(1:npti,:) / a_i_2d(1:npti,:) 
540         t_su_2d(1:npti,:)  =  zaTsfn(1:npti,:) / a_i_2d(1:npti,:) 
541      ELSEWHERE
542         h_i_2d(1:npti,:)  = 0._wp
543         t_su_2d(1:npti,:)  = rt0
544      END WHERE
545      !
546write(numout,*) "RSRH PE iceitd HH h_i", sum(h_i(:,:,:)); flush(numout)
547write(numout,*) "RSRH PE iceitd HH a_i", sum(a_i(:,:,:)); flush(numout)
548write(numout,*) "RSRH PE iceitd HH a_ip", sum(a_ip(:,:,:)); flush(numout)
549      CALL tab_2d_3d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i  )
550      CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i  )
551      CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i  )
552      CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s  )
553      CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i )
554      CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i )
555      CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
556      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
557      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
558write(numout,*) "RSRH PE iceitd II h_i", sum(h_i(:,:,:)); flush(numout)
559write(numout,*) "RSRH PE iceitd II a_i", sum(a_i(:,:,:)); flush(numout)
560write(numout,*) "RSRH PE iceitd II a_ip", sum(a_ip(:,:,:)); flush(numout)
561      !
562   END SUBROUTINE ice_itd_shiftice
563   
564
565   SUBROUTINE ice_itd_reb( kt )
566      !!------------------------------------------------------------------
567      !!                ***  ROUTINE ice_itd_reb ***
568      !!
569      !! ** Purpose : rebin - rebins thicknesses into defined categories
570      !!
571      !! ** Method  : If a category thickness is out of bounds, shift part (for down to top)
572      !!              or entire (for top to down) area, volume, and energy
573      !!              to the neighboring category
574      !!------------------------------------------------------------------
575      INTEGER , INTENT (in) ::   kt      ! Ocean time step
576      INTEGER ::   ji, jj, jl   ! dummy loop indices
577      !
578      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor           ! donor category index
579      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred
580      !!------------------------------------------------------------------
581      !
582      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 
583
584      jdonor(:,:) = 0
585      zdaice(:,:) = 0._wp
586      zdvice(:,:) = 0._wp
587      !
588      !                       !---------------------------------------
589      DO jl = 1, jpl-1        ! identify thicknesses that are too big
590         !                    !---------------------------------------
591         npti = 0   ;   nptidx(:) = 0
592         DO jj = 1, jpj
593            DO ji = 1, jpi
594               IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
595                  npti = npti + 1
596                  nptidx( npti ) = (jj - 1) * jpi + ji                 
597               ENDIF
598            END DO
599         END DO
600         !
601write(numout,*) "RSRH PE iceitd KK h_i",jl, sum(h_i(:,:,jl)); flush(numout)
602write(numout,*) "RSRH PE iceitd KK a_i",jl, sum(a_i(:,:,jl)); flush(numout)
603write(numout,*) "RSRH PE iceitd KK a_ip",jl, sum(a_ip(:,:,jl)); flush(numout)
604        !
605!!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)   )
606         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    )
607         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl)    )
608write(numout,*) "RSRH PE iceitd LL h_i",jl, sum(h_i(:,:,jl)); flush(numout)
609write(numout,*) "RSRH PE iceitd LL a_i",jl, sum(a_i(:,:,jl)); flush(numout)
610write(numout,*) "RSRH PE iceitd LL a_ip",jl, sum(a_ip(:,:jl)); flush(numout)
611         !
612         DO ji = 1, npti
613            jdonor(ji,jl)  = jl 
614            ! how much of a_i you send in cat sup is somewhat arbitrary
615!!clem: these do not work properly after a restart (I do not know why)
616!!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 
617!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
618!!clem: these do not work properly after a restart (I do not know why)
619!!            zdaice(ji,jl)  = a_i_1d(ji)
620!!            zdvice(ji,jl)  = v_i_1d(ji)
621!!clem: these are from UCL and work ok
622            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp
623            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
624         END DO
625         !
626         IF( npti > 0 ) THEN
627            CALL ice_itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1
628            ! Reset shift parameters
629            jdonor(1:npti,jl) = 0
630            zdaice(1:npti,jl) = 0._wp
631            zdvice(1:npti,jl) = 0._wp
632         ENDIF
633         !
634      END DO
635
636      !                       !-----------------------------------------
637      DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small
638write(numout,*) "RSRH PE iceitd LL for jl=",jl ;flush(numout)
639         !                    !-----------------------------------------
640         npti = 0 ; nptidx(:) = 0
641         DO jj = 1, jpj
642            DO ji = 1, jpi
643               IF( a_i(ji,jj,jl+1) > epsi10 .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN
644                  npti = npti + 1
645                  nptidx( npti ) = (jj - 1) * jpi + ji                 
646               ENDIF
647            END DO
648         END DO
649write(numout,*) "RSRH PE iceitd LL after 1st loop for jl=",jl ;flush(numout)
650write(numout,*) "RSRH PE iceitd get summing LL after 1st loop for jl=",jl ;flush(numout)
651         !
652write(numout,*) "RSRH PE iceitd MM pre h_i",jl, sum(h_i(:,:,jl)); flush(numout)
653write(numout,*) "RSRH PE iceitd MM pre a_i",jl, sum(a_i(:,:,jl)); flush(numout)
654write(numout,*) "RSRH PE iceitd MM pre a_ip",jl, sum(a_ip(:,:jl)); flush(numout)
655         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl+1)    ) ! jl+1 is ok
656         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl+1)    ) ! jl+1 is ok
657write(numout,*) "RSRH PE iceitd MM h_i",jl, sum(h_i(:,:,jl)); flush(numout)
658write(numout,*) "RSRH PE iceitd MM a_i",jl, sum(a_i(:,:,jl)); flush(numout)
659write(numout,*) "RSRH PE iceitd MM a_ip",jl, sum(a_ip(:,:jl)); flush(numout)
660         DO ji = 1, npti
661            jdonor(ji,jl) = jl + 1
662            zdaice(ji,jl) = a_i_1d(ji) 
663            zdvice(ji,jl) = v_i_1d(ji)
664         END DO
665write(numout,*) "RSRH PE iceitd MM after 2nd loop for jl=",jl ;flush(numout)
666         !
667         IF( npti > 0 ) THEN
668            CALL ice_itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl
669            ! Reset shift parameters
670            jdonor(1:npti,jl) = 0
671            zdaice(1:npti,jl) = 0._wp
672            zdvice(1:npti,jl) = 0._wp
673         ENDIF
674write(numout,*) "RSRH PE iceitd NN after 3rd loop for jl=",jl ;flush(numout)
675         !
676      END DO
677write(numout,*) "RSRH PE iceitd NN out of NN loop for jl=",jl ;flush(numout)
678      !
679   END SUBROUTINE ice_itd_reb
680
681   SUBROUTINE ice_itd_init
682      !!------------------------------------------------------------------
683      !!                ***  ROUTINE ice_itd_init ***
684      !!
685      !! ** Purpose :   Initializes the ice thickness distribution
686      !! ** Method  :   ...
687      !! ** input   :   Namelist namitd
688      !!-------------------------------------------------------------------
689      INTEGER  ::   jl            ! dummy loop index
690      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read
691      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
692      !
693      NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin
694      !!------------------------------------------------------------------
695      !
696      REWIND( numnam_ice_ref )      ! Namelist namitd in reference namelist : Parameters for ice
697      READ  ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901)
698901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namitd in reference namelist', lwp )
699
700      REWIND( numnam_ice_cfg )      ! Namelist namitd in configuration namelist : Parameters for ice
701      READ  ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 )
702902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namitd in configuration namelist', lwp )
703      IF(lwm) WRITE ( numoni, namitd )
704      !
705      IF(lwp) THEN                  ! control print
706         WRITE(numout,*)
707         WRITE(numout,*) 'ice_itd_init: Initialization of ice cat distribution '
708         WRITE(numout,*) '~~~~~~~~~~~~'
709         WRITE(numout,*) '   Namelist namitd: '
710         WRITE(numout,*) '      Ice categories are defined by a function of rn_himean**(-0.05)    ln_cat_hfn = ', ln_cat_hfn
711         WRITE(numout,*) '         mean ice thickness in the domain                               rn_himean  = ', rn_himean
712         WRITE(numout,*) '      Ice categories are defined by rn_catbnd                           ln_cat_usr = ', ln_cat_usr
713         WRITE(numout,*) '         ice categories boundaries (m)                                  rn_catbnd  = ', rn_catbnd 
714         WRITE(numout,*) '      minimum ice thickness                                             rn_himin   = ', rn_himin 
715      ENDIF
716      !
717      !-----------------------------------!
718      !  Thickness categories boundaries  !
719      !-----------------------------------!
720      !                             !== set the choice of ice categories ==!
721      ioptio = 0 
722      IF( ln_cat_hfn ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_cathfn    ;   ENDIF
723      IF( ln_cat_usr ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_catusr    ;   ENDIF
724      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_itd_init: choose one and only one ice categories boundaries' )
725      !
726      SELECT CASE( nice_catbnd )
727      !                                !------------------------!
728      CASE( np_cathfn )                ! h^(-alpha) function
729         !                             !------------------------!
730         zalpha = 0.05_wp
731         zhmax  = 3._wp * rn_himean
732         DO jl = 1, jpl
733            znum = jpl * ( zhmax+1 )**zalpha
734            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
735            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
736         END DO
737         !                             !------------------------!
738      CASE( np_catusr )                ! user defined
739         !                             !------------------------!
740         DO jl = 0, jpl
741            hi_max(jl) = rn_catbnd(jl)
742         END DO
743         !
744      END SELECT
745      !
746      DO jl = 1, jpl                ! mean thickness by category
747         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
748      END DO
749      !
750      hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
751      !
752      IF(lwp) WRITE(numout,*)
753      IF(lwp) WRITE(numout,*) '   ===>>>   resulting thickness category boundaries :'
754      IF(lwp) WRITE(numout,*) '            hi_max(:)= ', hi_max(0:jpl)
755      !
756   END SUBROUTINE ice_itd_init
757
758#else
759   !!----------------------------------------------------------------------
760   !!   Default option :         Empty module         NO ESIM sea-ice model
761   !!----------------------------------------------------------------------
762#endif
763
764   !!======================================================================
765END MODULE iceitd
Note: See TracBrowser for help on using the repository browser.