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.
limitd_th.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
  • Property svn:keywords set to Id
File size: 38.0 KB
Line 
1MODULE limitd_th
2   !!======================================================================
3   !!                       ***  MODULE limitd_th ***
4   !!   LIM3 ice model : ice thickness distribution: Thermodynamics
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' :                                   LIM3 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_itd_th_rem   :
16   !!   lim_itd_th_reb   :
17   !!   lim_itd_fitline  :
18   !!   lim_itd_shiftice :
19   !!----------------------------------------------------------------------
20   USE par_oce          ! ocean parameters
21   USE dom_oce          ! ocean domain
22   USE phycst           ! physical constants (ocean directory)
23   USE thd_ice          ! LIM-3 thermodynamic variables
24   USE ice              ! LIM-3 variables
25   USE limvar           ! LIM-3 variables
26   USE prtctl           ! Print control
27   USE in_out_manager   ! I/O manager
28   USE lib_mpp          ! MPP library
29   USE wrk_nemo         ! work arrays
30   USE lib_fortran      ! to use key_nosignedzero
31   USE limcons          ! conservation tests
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   lim_itd_th_rem
37   PUBLIC   lim_itd_th_reb
38
39   !!----------------------------------------------------------------------
40   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt )
47      !!------------------------------------------------------------------
48      !!                ***  ROUTINE lim_itd_th_rem ***
49      !!
50      !! ** Purpose :   computes the redistribution of ice thickness
51      !!              after thermodynamic growth of ice thickness
52      !!
53      !! ** Method  : Linear remapping
54      !!
55      !! References : W.H. Lipscomb, JGR 2001
56      !!------------------------------------------------------------------
57      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point
58      INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied
59      INTEGER , INTENT (in) ::   kt      ! Ocean time step
60      !
61      INTEGER  ::   ji, jj, jl     ! dummy loop index
62      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji
63      INTEGER  ::   nd             ! local integer
64      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars
65      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
66      REAL(wp) ::   zx3       
67      CHARACTER (len = 15) :: fieldid
68
69      INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor   ! donor category index
70
71      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdhice      ! ice thickness increment
72      REAL(wp), POINTER, DIMENSION(:,:,:) ::   g0          ! coefficients for fitting the line of the ITD
73      REAL(wp), POINTER, DIMENSION(:,:,:) ::   g1          ! coefficients for fitting the line of the ITD
74      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hL          ! left boundary for the ITD for each thickness
75      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hR          ! left boundary for the ITD for each thickness
76      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zht_i_b     ! old ice thickness
77      REAL(wp), POINTER, DIMENSION(:,:,:) ::   dummy_es
78      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice          ! local increment of ice area and volume
79      REAL(wp), POINTER, DIMENSION(:)     ::   zvetamin, zvetamax      ! maximum values for etas
80      INTEGER , POINTER, DIMENSION(:)     ::   nind_i, nind_j          ! compressed indices for i/j directions
81      INTEGER                             ::   nbrem                   ! number of cells with ice to transfer
82      REAL(wp)                            ::   zslope                  ! used to compute local thermodynamic "speeds"
83      REAL(wp), POINTER, DIMENSION(:,:)   ::   zhb0, zhb1              ! category boundaries for thinnes categories
84      REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_i_init, vt_i_final   !  ice volume summed over categories
85      REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_s_init, vt_s_final   !  snow volume summed over categories
86      REAL(wp), POINTER, DIMENSION(:,:)   ::   et_i_init, et_i_final   !  ice energy summed over categories
87      REAL(wp), POINTER, DIMENSION(:,:)   ::   et_s_init, et_s_final   !  snow energy summed over categories
88      INTEGER , POINTER, DIMENSION(:,:)   ::   zremap_flag      ! compute remapping or not ????
89      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhbnew           ! new boundaries of ice categories
90      !!------------------------------------------------------------------
91
92      CALL wrk_alloc( jpi,jpj, zremap_flag )
93      CALL wrk_alloc( jpi,jpj,jpl-1, zdonor )
94      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es )
95      CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )   
96      CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )   
97      CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )   
98      CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
99      CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final )
100
101      !!----------------------------------------------------------------------------------------------
102      !! 0) Conservation checkand changes in each ice category
103      !!----------------------------------------------------------------------------------------------
104      IF( con_i ) THEN
105         CALL lim_column_sum (jpl,   v_i, vt_i_init)
106         CALL lim_column_sum (jpl,   v_s, vt_s_init)
107         CALL lim_column_sum_energy (jpl, nlay_i,   e_i, et_i_init)
108         DO jl = 1, jpl
109!$OMP PARALLEL DO schedule(static) private(jj,ji)
110            DO jj = 1, jpj
111               DO ji = 1, jpi
112                  dummy_es(ji,jj,jl) = e_s(ji,jj,1,jl)
113               END DO
114            END DO
115         END DO
116         CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init)
117      ENDIF
118
119      !!----------------------------------------------------------------------------------------------
120      !! 1) Compute thickness and changes in each ice category
121      !!----------------------------------------------------------------------------------------------
122      IF( kt == nit000 .AND. lwp) THEN
123         WRITE(numout,*)
124         WRITE(numout,*) 'lim_itd_th_rem  : Remapping the ice thickness distribution'
125         WRITE(numout,*) '~~~~~~~~~~~~~~~'
126         WRITE(numout,*) ' klbnd :       ', klbnd
127         WRITE(numout,*) ' kubnd :       ', kubnd
128      ENDIF
129
130!$OMP PARALLEL
131      DO jl = 1, jpl
132!$OMP DO schedule(static) private(jj,ji)
133         DO jj = 1, jpj
134            DO ji = 1, jpi
135               zdhice(ji,jj,jl) = 0._wp
136            END DO
137         END DO
138      END DO
139      DO jl = klbnd, kubnd
140!$OMP DO schedule(static) private(jj,ji,rswitch)
141         DO jj = 1, jpj
142            DO ji = 1, jpi
143               rswitch           = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) )     !0 if no ice and 1 if yes
144               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch
145               rswitch           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) )
146               zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch
147               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) ! clem: useless IF statement?
148            END DO
149         END DO
150      END DO
151
152      !-----------------------------------------------------------------------------------------------
153      !  2) Compute fractional ice area in each grid cell
154      !-----------------------------------------------------------------------------------------------
155!$OMP DO schedule(static) private(jj,ji)
156      DO jj = 1, jpj
157         DO ji = 1, jpi
158            at_i(ji,jj) = 0._wp
159         END DO
160      END DO
161      DO jl = klbnd, kubnd
162!$OMP DO schedule(static) private(jj,ji)
163         DO jj = 1, jpj
164            DO ji = 1, jpi
165               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
166            END DO
167         END DO
168      END DO
169!$OMP END PARALLEL
170
171      !-----------------------------------------------------------------------------------------------
172      !  3) Identify grid cells with ice
173      !-----------------------------------------------------------------------------------------------
174      nbrem = 0
175      DO jj = 1, jpj
176         DO ji = 1, jpi
177            IF ( at_i(ji,jj) > epsi10 ) THEN
178               nbrem         = nbrem + 1
179               nind_i(nbrem) = ji
180               nind_j(nbrem) = jj
181               zremap_flag(ji,jj) = 1
182            ELSE
183               zremap_flag(ji,jj) = 0
184            ENDIF
185         END DO
186      END DO
187
188      !-----------------------------------------------------------------------------------------------
189      !  4) Compute new category boundaries
190      !-----------------------------------------------------------------------------------------------
191      !- 4.1 Compute category boundaries
192!$OMP PARALLEL
193      DO jl = 0, jpl
194!$OMP DO schedule(static) private(jj,ji)
195         DO jj = 1, jpj
196            DO ji = 1, jpi
197               zhbnew(ji,jj,jl) = 0._wp
198            END DO
199         END DO
200      END DO
201
202      DO jl = klbnd, kubnd - 1
203!$OMP DO schedule(static) private(ji,ii,ij,zslope)
204         DO ji = 1, nbrem
205            ii = nind_i(ji)
206            ij = nind_j(ji)
207            !
208            zhbnew(ii,ij,jl) = hi_max(jl)
209            IF    ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN
210               !interpolate between adjacent category growth rates
211               zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) )
212               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) )
213            ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN
214               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl)
215            ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN
216               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1)
217            ENDIF
218         END DO
219
220         !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness
221!$OMP DO schedule(static) private(ji,ii,ij)
222         DO ji = 1, nbrem
223            ii = nind_i(ji)
224            ij = nind_j(ji)
225
226            ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible
227            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
228            IF    ( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN
229               zremap_flag(ii,ij) = 0
230            ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN
231               zremap_flag(ii,ij) = 0
232            ENDIF
233
234            !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 
235            IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0
236            IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0
237            ! clem bug: why is not the following instead?
238            !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0
239            !!IF( zhbnew(ii,ij,jl) > hi_max(jl  ) ) zremap_flag(ii,ij) = 0
240 
241         END DO
242
243      END DO
244!$OMP END PARALLEL
245
246      !-----------------------------------------------------------------------------------------------
247      !  5) Identify cells where ITD is to be remapped
248      !-----------------------------------------------------------------------------------------------
249      nbrem = 0
250      DO jj = 1, jpj
251         DO ji = 1, jpi
252            IF( zremap_flag(ji,jj) == 1 ) THEN
253               nbrem         = nbrem + 1
254               nind_i(nbrem) = ji
255               nind_j(nbrem) = jj
256            ENDIF
257         END DO
258      END DO 
259
260      !-----------------------------------------------------------------------------------------------
261      !  6) Fill arrays with lowermost / uppermost boundaries of 'new' categories
262      !-----------------------------------------------------------------------------------------------
263!$OMP PARALLEL DO schedule(static) private(jj,ji)
264      DO jj = 1, jpj
265         DO ji = 1, jpi
266            zhb0(ji,jj) = hi_max(0)
267            zhb1(ji,jj) = hi_max(1)
268
269            IF( a_i(ji,jj,kubnd) > epsi10 ) THEN
270               zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) )
271            ELSE
272!clem bug               zhbnew(ji,jj,kubnd) = hi_max(kubnd) 
273               zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway
274            ENDIF
275
276            ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible
277            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
278            IF    ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) )  THEN
279               zremap_flag(ji,jj) = 0
280            ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) )  THEN
281               zremap_flag(ji,jj) = 0
282            ENDIF
283
284         END DO
285      END DO
286
287      !-----------------------------------------------------------------------------------------------
288      !  7) Compute g(h)
289      !-----------------------------------------------------------------------------------------------
290      !- 7.1 g(h) for category 1 at start of time step
291      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   &
292         &                  hR(:,:,klbnd), zremap_flag )
293
294      !- 7.2 Area lost due to melting of thin ice (first category,  klbnd)
295!$OMP PARALLEL DO schedule(static) private(ji,ii,ij,zdh0,zetamax,zx1,zx2,zda0,zdamax)
296      DO ji = 1, nbrem
297         ii = nind_i(ji) 
298         ij = nind_j(ji) 
299
300         IF( a_i(ii,ij,klbnd) > epsi10 ) THEN
301
302            zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category
303
304            IF( zdh0 < 0.0 ) THEN      !remove area from category 1
305               zdh0 = MIN( -zdh0, hi_max(klbnd) )
306               !Integrate g(1) from 0 to dh0 to estimate area melted
307               zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd)
308
309               IF( zetamax > 0.0 ) THEN
310                  zx1    = zetamax
311                  zx2    = 0.5 * zetamax * zetamax 
312                  zda0   = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1                        ! ice area removed
313                  zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i               
314                  zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting
315                                                                                                !     of thin ice (zdamax > 0)
316                  ! Remove area, conserving volume
317                  ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 )
318                  a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0
319                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ?
320               ENDIF
321
322            ELSE ! if ice accretion zdh0 > 0
323               ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
324               zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 
325            ENDIF
326
327         ENDIF
328
329      END DO
330
331      !- 7.3 g(h) for each thickness category 
332      DO jl = klbnd, kubnd
333         CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),  &
334            &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag )
335      END DO
336
337      !-----------------------------------------------------------------------------------------------
338      !  8) Compute area and volume to be shifted across each boundary
339      !-----------------------------------------------------------------------------------------------
340
341!$OMP PARALLEL
342      DO jl = klbnd, kubnd - 1
343!$OMP DO schedule(static) private(jj,ji)
344         DO jj = 1, jpj
345            DO ji = 1, jpi
346               zdonor(ji,jj,jl) = 0
347               zdaice(ji,jj,jl) = 0.0
348               zdvice(ji,jj,jl) = 0.0
349            END DO
350         END DO
351
352!$OMP DO schedule(static) private(ji,ii,ij,zetamax,zetamin,zx1,zwk1,zwk2,zx2,zx3,nd)
353         DO ji = 1, nbrem
354            ii = nind_i(ji)
355            ij = nind_j(ji)
356
357            IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1
358               ! left and right integration limits in eta space
359               zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)
360               zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl)
361               zdonor(ii,ij,jl) = jl
362
363            ELSE                                    ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl
364               ! left and right integration limits in eta space
365               zvetamin(ji) = 0.0
366               zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)
367               zdonor(ii,ij,jl) = jl + 1
368
369            ENDIF
370
371            zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin
372            zetamin = zvetamin(ji)
373
374            zx1  = zetamax - zetamin
375            zwk1 = zetamin * zetamin
376            zwk2 = zetamax * zetamax
377            zx2  = 0.5 * ( zwk2 - zwk1 )
378            zwk1 = zwk1 * zetamin
379            zwk2 = zwk2 * zetamax
380            zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 )
381            nd   = zdonor(ii,ij,jl)
382            zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1
383            zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd)
384
385         END DO
386      END DO
387!$OMP END PARALLEL
388
389      !!----------------------------------------------------------------------------------------------
390      !! 9) Shift ice between categories
391      !!----------------------------------------------------------------------------------------------
392      CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice )
393
394      !!----------------------------------------------------------------------------------------------
395      !! 10) Make sure ht_i >= minimum ice thickness hi_min
396      !!----------------------------------------------------------------------------------------------
397
398!$OMP PARALLEL DO schedule(static) private(ji,ii,ij)
399      DO ji = 1, nbrem
400         ii = nind_i(ji)
401         ij = nind_j(ji)
402         IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN
403            a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 
404            ht_i(ii,ij,1) = rn_himin
405         ENDIF
406      END DO
407
408      !!----------------------------------------------------------------------------------------------
409      !! 11) Conservation check
410      !!----------------------------------------------------------------------------------------------
411      IF ( con_i ) THEN
412         CALL lim_column_sum (jpl,   v_i, vt_i_final)
413         fieldid = ' v_i : limitd_th '
414         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
415
416         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, et_i_final)
417         fieldid = ' e_i : limitd_th '
418         CALL lim_cons_check (et_i_init, et_i_final, 1.0e-3, fieldid) 
419
420         CALL lim_column_sum (jpl,   v_s, vt_s_final)
421         fieldid = ' v_s : limitd_th '
422         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
423
424         DO jl = 1, jpl
425!$OMP PARALLEL DO schedule(static) private(jj,ji)
426            DO jj = 1, jpj
427               DO ji = 1, jpi
428                  dummy_es(ji,jj,jl) = e_s(ji,jj,1,jl)
429               END DO
430            END DO
431         END DO
432         CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_final)
433         fieldid = ' e_s : limitd_th '
434         CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid) 
435      ENDIF
436
437      CALL wrk_dealloc( jpi,jpj, zremap_flag )
438      CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor )
439      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es )
440      CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )   
441      CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )   
442      CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )   
443      CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
444      CALL wrk_dealloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final )
445
446   END SUBROUTINE lim_itd_th_rem
447
448
449   SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )
450      !!------------------------------------------------------------------
451      !!                ***  ROUTINE lim_itd_fitline ***
452      !!
453      !! ** Purpose :   fit g(h) with a line using area, volume constraints
454      !!
455      !! ** Method  :   Fit g(h) with a line, satisfying area and volume constraints.
456      !!              To reduce roundoff errors caused by large values of g0 and g1,
457      !!              we actually compute g(eta), where eta = h - hL, and hL is the
458      !!              left boundary.
459      !!------------------------------------------------------------------
460      INTEGER                     , INTENT(in   ) ::   num_cat      ! category index
461      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   HbL, HbR     ! left and right category boundaries
462      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   hice         ! ice thickness
463      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   g0, g1       ! coefficients in linear equation for g(eta)
464      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   hL           ! min value of range over which g(h) > 0
465      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   hR           ! max value of range over which g(h) > 0
466      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) ::   zremap_flag  !
467      !
468      INTEGER  ::   ji,jj        ! horizontal indices
469      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL)
470      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL)
471      REAL(wp) ::   zdhr         ! 1 / (hR - hL)
472      REAL(wp) ::   zwk1, zwk2   ! temporary variables
473      !!------------------------------------------------------------------
474      !
475!$OMP PARALLEL DO schedule(static) private(jj,ji,zh13,zh23,zdhr,zwk1,zwk2)
476      DO jj = 1, jpj
477         DO ji = 1, jpi
478            !
479            IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10   &
480               &                        .AND. hice(ji,jj)        > 0._wp )  THEN
481
482               ! Initialize hL and hR
483               hL(ji,jj) = HbL(ji,jj)
484               hR(ji,jj) = HbR(ji,jj)
485
486               ! Change hL or hR if hice falls outside central third of range
487               zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) )
488               zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) )
489
490               IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj)
491               ELSEIF( hice(ji,jj) > zh23 ) THEN   ;   hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj)
492               ENDIF
493
494               ! Compute coefficients of g(eta) = g0 + g1*eta
495               zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj))
496               zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr
497               zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr
498               g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 )
499               g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 )
500               !
501            ELSE  ! remap_flag = .false. or a_i < epsi10
502               hL(ji,jj) = 0._wp
503               hR(ji,jj) = 0._wp
504               g0(ji,jj) = 0._wp
505               g1(ji,jj) = 0._wp
506            ENDIF
507            !
508         END DO
509      END DO
510      !
511   END SUBROUTINE lim_itd_fitline
512
513
514   SUBROUTINE lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
515      !!------------------------------------------------------------------
516      !!                ***  ROUTINE lim_itd_shiftice ***
517      !!
518      !! ** Purpose :   shift ice across category boundaries, conserving everything
519      !!              ( area, volume, energy, age*vol, and mass of salt )
520      !!
521      !! ** Method  :
522      !!------------------------------------------------------------------
523      INTEGER                           , INTENT(in   ) ::   klbnd    ! Start thickness category index point
524      INTEGER                           , INTENT(in   ) ::   kubnd    ! End point on which the  the computation is applied
525      INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in   ) ::   zdonor   ! donor category index
526      REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdaice   ! ice area transferred across boundary
527      REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdvice   ! ice volume transferred across boundary
528
529      INTEGER ::   ji, jj, jl, jl2, jl1, jk   ! dummy loop indices
530      INTEGER ::   ii, ij                     ! indices when changing from 2D-1D is done
531
532      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zaTsfn
533      REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka            ! temporary array used here
534
535      REAL(wp) ::   zdvsnow, zdesnow   ! snow volume and energy transferred
536      REAL(wp) ::   zdeice             ! ice energy transferred
537      REAL(wp) ::   zdsm_vice          ! ice salinity times volume transferred
538      REAL(wp) ::   zdo_aice           ! ice age times volume transferred
539      REAL(wp) ::   zdaTsf             ! aicen*Tsfcn transferred
540
541      INTEGER, POINTER, DIMENSION(:) ::   nind_i, nind_j   ! compressed indices for i/j directions
542
543      INTEGER  ::   nbrem             ! number of cells with ice to transfer
544      !!------------------------------------------------------------------
545
546      CALL wrk_alloc( jpi,jpj,jpl, zaTsfn )
547      CALL wrk_alloc( jpi,jpj, zworka )
548      CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )
549
550      !----------------------------------------------------------------------------------------------
551      ! 1) Define a variable equal to a_i*T_su
552      !----------------------------------------------------------------------------------------------
553
554      DO jl = klbnd, kubnd
555!$OMP PARALLEL DO schedule(static) private(jj,ji)
556         DO jj = 1, jpj
557            DO ji = 1, jpi
558               zaTsfn(ji,jj,jl) = a_i(ji,jj,jl) * t_su(ji,jj,jl)
559            END DO
560         END DO
561      END DO
562
563      !-------------------------------------------------------------------------------
564      ! 2) Transfer volume and energy between categories
565      !-------------------------------------------------------------------------------
566
567      DO jl = klbnd, kubnd - 1
568         nbrem = 0
569         DO jj = 1, jpj
570            DO ji = 1, jpi
571               IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny
572                  nbrem = nbrem + 1
573                  nind_i(nbrem) = ji
574                  nind_j(nbrem) = jj
575               ENDIF
576            END DO
577         END DO
578
579!$OMP PARALLEL DO schedule(static) private(ji,ii,ij,jl1,jl2,rswitch,zdvsnow,zdesnow,zdo_aice,zdsm_vice,zdaTsf)
580         DO ji = 1, nbrem 
581            ii = nind_i(ji)
582            ij = nind_j(ji)
583
584            jl1 = zdonor(ii,ij,jl)
585            rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi10 ) )
586            zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch
587            IF( jl1 == jl) THEN   ;   jl2 = jl1+1
588            ELSE                  ;   jl2 = jl 
589            ENDIF
590
591            !--------------
592            ! Ice areas
593            !--------------
594            a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl)
595            a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl)
596
597            !--------------
598            ! Ice volumes
599            !--------------
600            v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 
601            v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl)
602
603            !--------------
604            ! Snow volumes
605            !--------------
606            zdvsnow        = v_s(ii,ij,jl1) * zworka(ii,ij)
607            v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow
608            v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow 
609
610            !--------------------
611            ! Snow heat content 
612            !--------------------
613            zdesnow            = e_s(ii,ij,1,jl1) * zworka(ii,ij)
614            e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow
615            e_s(ii,ij,1,jl2)   = e_s(ii,ij,1,jl2) + zdesnow
616
617            !--------------
618            ! Ice age
619            !--------------
620            zdo_aice           = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl)
621            oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice
622            oa_i(ii,ij,jl2)    = oa_i(ii,ij,jl2) + zdo_aice
623
624            !--------------
625            ! Ice salinity
626            !--------------
627            zdsm_vice          = smv_i(ii,ij,jl1) * zworka(ii,ij)
628            smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice
629            smv_i(ii,ij,jl2)   = smv_i(ii,ij,jl2) + zdsm_vice
630
631            !---------------------
632            ! Surface temperature
633            !---------------------
634            zdaTsf             = t_su(ii,ij,jl1) * zdaice(ii,ij,jl)
635            zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf
636            zaTsfn(ii,ij,jl2)  = zaTsfn(ii,ij,jl2) + zdaTsf 
637
638         END DO
639
640         !------------------
641         ! Ice heat content
642         !------------------
643
644         DO jk = 1, nlay_i
645!$OMP PARALLEL DO schedule(static) private(ji,ii,ij,jl1,jl2,zdeice)
646            DO ji = 1, nbrem
647               ii = nind_i(ji)
648               ij = nind_j(ji)
649
650               jl1 = zdonor(ii,ij,jl)
651               IF (jl1 == jl) THEN
652                  jl2 = jl+1
653               ELSE             ! n1 = n+1
654                  jl2 = jl 
655               ENDIF
656
657               zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij)
658               e_i(ii,ij,jk,jl1) =  e_i(ii,ij,jk,jl1) - zdeice
659               e_i(ii,ij,jk,jl2) =  e_i(ii,ij,jk,jl2) + zdeice 
660            END DO
661         END DO
662
663      END DO                   ! boundaries, 1 to ncat-1
664
665      !-----------------------------------------------------------------
666      ! Update ice thickness and temperature
667      !-----------------------------------------------------------------
668
669      DO jl = klbnd, kubnd
670!$OMP PARALLEL DO schedule(static) private(jj,ji)
671         DO jj = 1, jpj
672            DO ji = 1, jpi 
673               IF ( a_i(ji,jj,jl) > epsi10 ) THEN
674                  ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl) 
675                  t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 
676               ELSE
677                  ht_i(ji,jj,jl)  = 0._wp
678                  t_su(ji,jj,jl)  = rt0
679               ENDIF
680            END DO
681         END DO
682      END DO
683      !
684      CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn )
685      CALL wrk_dealloc( jpi,jpj, zworka )
686      CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )
687      !
688   END SUBROUTINE lim_itd_shiftice
689   
690
691   SUBROUTINE lim_itd_th_reb( klbnd, kubnd )
692      !!------------------------------------------------------------------
693      !!                ***  ROUTINE lim_itd_th_reb ***
694      !!
695      !! ** Purpose : rebin - rebins thicknesses into defined categories
696      !!
697      !! ** Method  :
698      !!------------------------------------------------------------------
699      INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point
700      INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied
701      !
702      INTEGER ::   ji,jj, jl   ! dummy loop indices
703      INTEGER ::   zshiftflag          ! = .true. if ice must be shifted
704      CHARACTER (len = 15) :: fieldid
705
706      INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor           ! donor category index
707      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice   ! ice area and volume transferred
708
709      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final   ! ice volume summed over categories
710      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   ! snow volume summed over categories
711      !!------------------------------------------------------------------
712     
713      CALL wrk_alloc( jpi,jpj,jpl, zdonor )   ! interger
714      CALL wrk_alloc( jpi,jpj,jpl, zdaice, zdvice )
715      CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final )
716      !     
717      IF( con_i ) THEN                 ! conservation check
718         CALL lim_column_sum (jpl,   v_i, vt_i_init)
719         CALL lim_column_sum (jpl,   v_s, vt_s_init)
720      ENDIF
721
722      !
723      !------------------------------------------------------------------------------
724      ! 1) Compute ice thickness.
725      !------------------------------------------------------------------------------
726!$OMP PARALLEL
727      DO jl = klbnd, kubnd
728!$OMP DO schedule(static) private(jj,ji,rswitch)
729         DO jj = 1, jpj
730            DO ji = 1, jpi 
731               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
732               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
733            END DO
734         END DO
735      END DO
736
737      !------------------------------------------------------------------------------
738      ! 2) If a category thickness is not in bounds, shift the
739      ! entire area, volume, and energy to the neighboring category
740      !------------------------------------------------------------------------------
741      !-------------------------
742      ! Initialize shift arrays
743      !-------------------------
744      DO jl = klbnd, kubnd
745!$OMP DO schedule(static) private(jj,ji)
746         DO jj = 1, jpj
747            DO ji = 1, jpi
748               zdonor(ji,jj,jl) = 0
749               zdaice(ji,jj,jl) = 0._wp
750               zdvice(ji,jj,jl) = 0._wp
751            END DO
752         END DO
753      END DO
754!$OMP END PARALLEL
755
756      !-------------------------
757      ! Move thin categories up
758      !-------------------------
759
760      DO jl = klbnd, kubnd - 1  ! loop over category boundaries
761
762         !---------------------------------------
763         ! identify thicknesses that are too big
764         !---------------------------------------
765         zshiftflag = 0
766
767!$OMP PARALLEL DO schedule(static) private(jj,ji) REDUCTION(MAX:zshiftflag)
768         DO jj = 1, jpj 
769            DO ji = 1, jpi 
770               IF( a_i(ji,jj,jl) > epsi10 .AND. ht_i(ji,jj,jl) > hi_max(jl) ) THEN
771                  zshiftflag        = 1
772                  zdonor(ji,jj,jl)  = jl 
773                  ! begin TECLIM change
774                  !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * 0.5_wp
775                  !zdvice(ji,jj,jl)  = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp
776                  ! end TECLIM change
777                  ! clem: how much of a_i you send in cat sup is somewhat arbitrary
778                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl) 
779                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 )
780               ENDIF
781            END DO
782         END DO
783         IF(lk_mpp)   CALL mpp_max( zshiftflag )
784
785         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
786            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
787            ! Reset shift parameters
788!$OMP PARALLEL DO schedule(static) private(jj,ji)
789            DO jj = 1, jpj
790               DO ji = 1, jpi
791                  zdonor(ji,jj,jl) = 0
792                  zdaice(ji,jj,jl) = 0._wp
793                  zdvice(ji,jj,jl) = 0._wp
794               END DO
795            END DO
796         ENDIF
797         !
798      END DO
799
800      !----------------------------
801      ! Move thick categories down
802      !----------------------------
803
804      DO jl = kubnd - 1, 1, -1       ! loop over category boundaries
805
806         !-----------------------------------------
807         ! Identify thicknesses that are too small
808         !-----------------------------------------
809         zshiftflag = 0
810
811!$OMP PARALLEL DO schedule(static) private(jj,ji) REDUCTION(MAX:zshiftflag)
812         DO jj = 1, jpj
813            DO ji = 1, jpi
814               IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN
815                  zshiftflag = 1
816                  zdonor(ji,jj,jl) = jl + 1
817                  zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 
818                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)
819               ENDIF
820            END DO
821         END DO
822
823         IF(lk_mpp)   CALL mpp_max( zshiftflag )
824         
825         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories
826            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
827            ! Reset shift parameters
828!$OMP PARALLEL DO schedule(static) private(jj,ji)
829            DO jj = 1, jpj
830               DO ji = 1, jpi
831                  zdonor(ji,jj,jl) = 0
832                  zdaice(ji,jj,jl) = 0._wp
833                  zdvice(ji,jj,jl) = 0._wp
834               END DO
835            END DO
836         ENDIF
837
838      END DO
839
840      !------------------------------------------------------------------------------
841      ! 3) Conservation check
842      !------------------------------------------------------------------------------
843
844      IF( con_i ) THEN
845         CALL lim_column_sum (jpl,   v_i, vt_i_final)
846         fieldid = ' v_i : limitd_reb '
847         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
848
849         CALL lim_column_sum (jpl,   v_s, vt_s_final)
850         fieldid = ' v_s : limitd_reb '
851         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
852      ENDIF
853      !
854      CALL wrk_dealloc( jpi,jpj,jpl, zdonor )
855      CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice )
856      CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final )
857
858   END SUBROUTINE lim_itd_th_reb
859
860#else
861   !!----------------------------------------------------------------------
862   !!   Default option            Dummy module         NO LIM sea-ice model
863   !!----------------------------------------------------------------------
864CONTAINS
865   SUBROUTINE lim_itd_th_rem
866   END SUBROUTINE lim_itd_th_rem
867   SUBROUTINE lim_itd_fitline
868   END SUBROUTINE lim_itd_fitline
869   SUBROUTINE lim_itd_shiftice
870   END SUBROUTINE lim_itd_shiftice
871   SUBROUTINE lim_itd_th_reb
872   END SUBROUTINE lim_itd_th_reb
873#endif
874   !!======================================================================
875END MODULE limitd_th
Note: See TracBrowser for help on using the repository browser.