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 branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 9 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

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