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.
limtrp.F90 in branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 5053

Last change on this file since 5053 was 5053, checked in by clem, 9 years ago

LIM3 cleaning continues. No change in the physics besides the introduction of the monocategory sea ice
  • Property svn:keywords set to Id
File size: 28.7 KB
Line 
1MODULE limtrp
2   !!======================================================================
3   !!                       ***  MODULE limtrp   ***
4   !! LIM transport ice model : sea-ice advection/diffusion
5   !!======================================================================
6   !! History : LIM-2 ! 2000-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet)  Original code
7   !!            3.0  ! 2005-11 (M. Vancoppenolle)   Multi-layer sea ice, salinity variations
8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                      LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_trp      : advection/diffusion process of sea ice
15   !!----------------------------------------------------------------------
16   USE phycst         ! physical constant
17   USE dom_oce        ! ocean domain
18   USE sbc_oce        ! ocean surface boundary condition
19   USE par_ice        ! ice parameter
20   USE dom_ice        ! ice domain
21   USE ice            ! ice variables
22   USE limadv         ! ice advection
23   USE limhdf         ! ice horizontal diffusion
24   USE limvar         !
25   !
26   USE in_out_manager ! I/O manager
27   USE lbclnk         ! lateral boundary conditions -- MPP exchanges
28   USE lib_mpp        ! MPP library
29   USE wrk_nemo       ! work arrays
30   USE prtctl         ! Print control
31   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
32   USE timing         ! Timing
33   USE limcons        ! conservation tests
34   USE limctl         ! control prints
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   lim_trp    ! called by sbcice_lim
40
41   INTEGER  ::   ncfl                 ! number of ice time step with CFL>1/2 
42
43   !! * Substitution
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE lim_trp( kt ) 
53      !!-------------------------------------------------------------------
54      !!                   ***  ROUTINE lim_trp ***
55      !!                   
56      !! ** purpose : advection/diffusion process of sea ice
57      !!
58      !! ** method  : variables included in the process are scalar,   
59      !!     other values are considered as second order.
60      !!     For advection, a second order Prather scheme is used. 
61      !!
62      !! ** action :
63      !!---------------------------------------------------------------------
64      INTEGER, INTENT(in) ::   kt   ! number of iteration
65      !
66      INTEGER  ::   ji, jj, jk, jl, jt   ! dummy loop indices
67      INTEGER  ::   initad                  ! number of sub-timestep for the advection
68      REAL(wp) ::   zcfl , zusnit           !   -      -
69      CHARACTER(len=80) :: cltmp
70      !
71      REAL(wp), POINTER, DIMENSION(:,:)      ::   zsm, zs0at
72      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e
73      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ow
74      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e
75      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold           ! old ice volume...
76      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness
77      REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold   ! old concentration, enthalpies
78      REAL(wp) :: zdv, zvi, zvs, zsmv, zes, zei
79      !
80      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
81      !!---------------------------------------------------------------------
82      IF( nn_timing == 1 )  CALL timing_start('limtrp')
83
84      CALL wrk_alloc( jpi,jpj,           zsm, zs0at, zatold, zeiold, zesold )
85      CALL wrk_alloc( jpi,jpj,jpl,       zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e )
86      CALL wrk_alloc( jpi,jpj,1,         zs0ow )
87      CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, zs0e )
88      CALL wrk_alloc( jpi,jpj,jpl,       zhimax, zviold, zvsold )
89
90      IF( numit == nstart .AND. lwp ) THEN
91         WRITE(numout,*)
92         IF( ln_limdyn ) THEN   ;   WRITE(numout,*) 'lim_trp : Ice transport '
93         ELSE                   ;   WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
94         ENDIF
95         WRITE(numout,*) '~~~~~~~~~~~~'
96         ncfl = 0                ! nb of time step with CFL > 1/2
97      ENDIF
98
99      zsm(:,:) = area(:,:)
100     
101      !                             !-------------------------------------!
102      IF( ln_limdyn ) THEN          !   Advection of sea ice properties   !
103         !                          !-------------------------------------!
104
105         ! conservation test
106         IF( ln_limdiahsb )   CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
107
108         ! mass and salt flux init
109         zviold(:,:,:) = v_i(:,:,:)
110         zeiold(:,:)   = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 
111         zesold(:,:)   = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 
112
113         !--- Thickness correction init. -------------------------------
114         CALL lim_var_glo2eqv
115         zatold(:,:) = SUM( a_i(:,:,:), dim=3 )
116         !---------------------------------------------------------------------
117         ! Record max of the surrounding ice thicknesses for correction in limupdate
118         ! in case advection creates ice too thick.
119         !---------------------------------------------------------------------
120         zhimax(:,:,:) = ht_i(:,:,:)
121         DO jl = 1, jpl
122            DO jj = 2, jpjm1
123               DO ji = 2, jpim1
124                  zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) )
125                  !zhimax(ji,jj,jl) = ( ht_i(ji  ,jj  ,jl) * tmask(ji,  jj  ,1) + ht_i(ji-1,jj-1,jl) * tmask(ji-1,jj-1,1) + ht_i(ji+1,jj+1,jl) * tmask(ji+1,jj+1,1) &
126                  !     &             + ht_i(ji-1,jj  ,jl) * tmask(ji-1,jj  ,1) + ht_i(ji  ,jj-1,jl) * tmask(ji  ,jj-1,1) &
127                  !     &             + ht_i(ji+1,jj  ,jl) * tmask(ji+1,jj  ,1) + ht_i(ji  ,jj+1,jl) * tmask(ji  ,jj+1,1) &
128                  !     &             + ht_i(ji-1,jj+1,jl) * tmask(ji-1,jj+1,1) + ht_i(ji+1,jj-1,jl) * tmask(ji+1,jj-1,1) )
129               END DO
130            END DO
131            CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
132         END DO
133         
134         !=============================!
135         !==      Prather scheme     ==!
136         !=============================!
137
138         ! If ice drift field is too fast, use an appropriate time step for advection.         
139         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) )         ! CFL test for stability
140         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) )
141         IF(lk_mpp )   CALL mpp_max( zcfl )
142!!gm more readability:
143!         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
144!         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
145!         ENDIF
146!!gm end
147         initad = 1 + NINT( MAX( 0._wp, SIGN( 1._wp, zcfl-0.5 ) ) )
148         zusnit = 1.0 / REAL( initad ) 
149         IF( zcfl > 0.5_wp .AND. lwp )   ncfl = ncfl + 1
150         IF( numit == nlast .AND. lwp ) THEN
151            IF( ncfl > 0 ) THEN   
152             WRITE(cltmp,'(i6.1)') ncfl
153             CALL ctl_stop('STOP',TRIM(cltmp) )
154               CALL ctl_warn( 'lim_trp: ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
155            ELSE
156               WRITE(numout,*) 'lim_trp : CFL criteria for ice advection is always smaller than 1/2 '
157            ENDIF
158         ENDIF
159
160         !-------------------------
161         ! transported fields                                       
162         !-------------------------
163         zs0ow(:,:,1) = ato_i(:,:) * area(:,:)               ! Open water area
164         DO jl = 1, jpl
165            zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume
166            zs0ice(:,:,jl)   = v_i  (:,:,jl) * area(:,:)    ! Ice  volume
167            zs0a  (:,:,jl)   = a_i  (:,:,jl) * area(:,:)    ! Ice area
168            zs0sm (:,:,jl)   = smv_i(:,:,jl) * area(:,:)    ! Salt content
169            zs0oi (:,:,jl)   = oa_i (:,:,jl) * area(:,:)    ! Age content
170            zs0c0 (:,:,jl)   = e_s  (:,:,1,jl)              ! Snow heat content
171            zs0e  (:,:,:,jl) = e_i  (:,:,:,jl)              ! Ice  heat content
172         END DO
173
174
175         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
176            DO jt = 1, initad
177               CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &             !--- ice open water area
178                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
179               CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &
180                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
181               DO jl = 1, jpl
182                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
183                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
184                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
185                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
186                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
187                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
188                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
189                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
190                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
191                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
192                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
193                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
194                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
195                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
196                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &
197                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
198                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
199                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
200                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
201                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
202                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
203                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
204                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
205                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
206                  DO jk = 1, nlay_i                                                                !--- ice heat contents ---
207                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
208                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
209                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
210                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
211                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
212                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
213                  END DO
214               END DO
215            END DO
216         ELSE
217            DO jt = 1, initad
218               CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &             !--- ice open water area
219                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
220               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &
221                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
222               DO jl = 1, jpl
223                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
224                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
225                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
226                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
227                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
228                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
229                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
230                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
231                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
232                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
233                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
234                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
235
236                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
237                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
238                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &
239                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
240                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
241                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
242                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &
243                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
244                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
245                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
246                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
247                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
248                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
249                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
250                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
251                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
252                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
253                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
254                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
255                  END DO
256               END DO
257            END DO
258         ENDIF
259
260         !-------------------------------------------
261         ! Recover the properties from their contents
262         !-------------------------------------------
263         ato_i(:,:) = zs0ow(:,:,1) / area(:,:)
264         DO jl = 1, jpl
265            v_i  (:,:,jl)   = zs0ice(:,:,jl) / area(:,:)
266            v_s  (:,:,jl)   = zs0sn (:,:,jl) / area(:,:)
267            smv_i(:,:,jl)   = zs0sm (:,:,jl) / area(:,:)
268            oa_i (:,:,jl)   = zs0oi (:,:,jl) / area(:,:)
269            a_i  (:,:,jl)   = zs0a  (:,:,jl) / area(:,:)
270            e_s  (:,:,1,jl) = zs0c0 (:,:,jl)             
271            e_i  (:,:,:,jl) = zs0e  (:,:,:,jl)           
272         END DO
273
274         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
275         DO jl = 2, jpl
276            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
277         END DO
278
279         !------------------------------------------------------------------------------!
280         ! 4) Diffusion of Ice fields                 
281         !------------------------------------------------------------------------------!
282
283         !
284         !--------------------------------
285         !  diffusion of open water area
286         !--------------------------------
287         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
288         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
289            DO ji = 1 , fs_jpim1   ! vector opt.
290               pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
291                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
292               pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
293                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
294            END DO
295         END DO
296         !
297         CALL lim_hdf( ato_i (:,:) )   ! Diffusion
298
299         !------------------------------------
300         !  Diffusion of other ice variables
301         !------------------------------------
302         DO jl = 1, jpl
303         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
304            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
305               DO ji = 1 , fs_jpim1   ! vector opt.
306                  pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,jl) ) ) )   &
307                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
308                  pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj  ,jl) ) ) )   &
309                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
310               END DO
311            END DO
312
313            CALL lim_hdf( v_i  (:,:,  jl) )
314            CALL lim_hdf( v_s  (:,:,  jl) )
315            CALL lim_hdf( smv_i(:,:,  jl) )
316            CALL lim_hdf( oa_i (:,:,  jl) )
317            CALL lim_hdf( a_i  (:,:,  jl) )
318            CALL lim_hdf( e_s  (:,:,1,jl) )
319            DO jk = 1, nlay_i
320               CALL lim_hdf( e_i(:,:,jk,jl) )
321            END DO
322         END DO
323
324         !------------------------------------------------------------------------------!
325         ! 5) Update and limit ice properties after transport                           
326         !------------------------------------------------------------------------------!
327
328!!gm & cr   :  MAX should not be active if adv scheme is positive !
329         !--------------------------------------------------
330         ! 5.1) Recover mean values over the grid squares.
331         !--------------------------------------------------
332         DO jl = 1, jpl
333            DO jj = 1, jpj
334               DO ji = 1, jpi
335                  v_s  (ji,jj,jl)   = MAX( 0._wp, v_s  (ji,jj,jl) )
336                  v_i  (ji,jj,jl)   = MAX( 0._wp, v_i  (ji,jj,jl) )
337                  smv_i(ji,jj,jl)   = MAX( 0._wp, smv_i(ji,jj,jl) )
338                  oa_i (ji,jj,jl)   = MAX( 0._wp, oa_i (ji,jj,jl) )
339                  a_i  (ji,jj,jl)   = MAX( 0._wp, a_i  (ji,jj,jl) )
340                  e_s  (ji,jj,1,jl) = MAX( 0._wp, e_s  (ji,jj,1,jl) )
341               END DO
342            END DO
343         END DO
344         DO jl = 1, jpl
345            DO jk = 1, nlay_i
346               DO jj = 1, jpj
347                  DO ji = 1, jpi
348                     e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) )
349                  END DO
350               END DO
351            END DO
352         END DO
353!!gm & cr
354
355         ! zap small areas
356         CALL lim_var_zapsmall
357
358         !--- Thickness correction in case too high --------------------------------------------------------
359         CALL lim_var_glo2eqv
360         DO jl = 1, jpl
361            DO jj = 1, jpj
362               DO ji = 1, jpi
363
364                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
365                     zvi  = v_i  (ji,jj,jl)
366                     zvs  = v_s  (ji,jj,jl)
367                     zsmv = smv_i(ji,jj,jl)
368                     zes  = e_s  (ji,jj,1,jl)
369                     zei  = SUM( e_i(ji,jj,1:nlay_i,jl) )
370                     zdv  = v_i(ji,jj,jl) - zviold(ji,jj,jl)   
371                     
372                     rswitch = 1._wp
373                     IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
374                        & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                         
375                        ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) )
376                        rswitch        = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )
377                        a_i(ji,jj,jl)  = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )
378                     ELSE
379                        ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) )
380                        rswitch        = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )
381                        a_i(ji,jj,jl)  = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )
382                     ENDIF
383
384                     ! small correction due to *rswitch for a_i
385                     v_i  (ji,jj,jl) = rswitch * v_i  (ji,jj,jl)
386                     v_s  (ji,jj,jl) = rswitch * v_s  (ji,jj,jl)
387                     smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl)
388                     e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl)
389                     e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
390
391                     ! Update mass fluxes
392                     wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
393                     wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
394                     sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
395                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
396                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
397                  ENDIF
398
399               END DO
400            END DO
401         END DO
402         ! -------------------------------------------------
403         
404         !--------------------------------------
405         ! Impose a_i < amax in mono-category
406         !--------------------------------------
407         !
408         IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2
409            DO jj = 1, jpj
410               DO ji = 1, jpi
411                  a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), amax )
412               END DO
413            END DO
414         ENDIF
415
416         ! --- diags ---
417         DO jj = 1, jpj
418            DO ji = 1, jpi
419               diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice
420               diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice
421
422               diag_trp_vi(ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice
423               diag_trp_vs(ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice
424            END DO
425         END DO
426
427         ! --- agglomerate variables -----------------
428         vt_i (:,:) = 0._wp
429         vt_s (:,:) = 0._wp
430         at_i (:,:) = 0._wp
431         !
432         DO jl = 1, jpl
433            DO jj = 1, jpj
434               DO ji = 1, jpi
435                  !
436                  vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
437                  vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
438                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
439               END DO
440            END DO
441         END DO
442         ! -------------------------------------------------
443
444         ! open water
445         DO jj = 1, jpj
446            DO ji = 1, jpi
447               ! open water = 1 if at_i=0
448               rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
449               ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
450            END DO
451         END DO     
452
453         ! conservation test
454         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
455
456      ENDIF
457
458      CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting
459
460      ! -------------------------------------------------
461      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print
462      ! -------------------------------------------------
463      IF(ln_ctl) THEN   ! Control print
464         CALL prt_ctl_info(' ')
465         CALL prt_ctl_info(' - Cell values : ')
466         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
467         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp  : cell area :')
468         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp  : at_i      :')
469         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp  : vt_i      :')
470         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp  : vt_s      :')
471         DO jl = 1, jpl
472            CALL prt_ctl_info(' ')
473            CALL prt_ctl_info(' - Category : ', ivar1=jl)
474            CALL prt_ctl_info('   ~~~~~~~~~~')
475            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_trp  : a_i      : ')
476            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_trp  : ht_i     : ')
477            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_trp  : ht_s     : ')
478            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_trp  : v_i      : ')
479            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_trp  : v_s      : ')
480            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_trp  : e_s      : ')
481            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_trp  : t_su     : ')
482            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_trp  : t_snow   : ')
483            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_trp  : sm_i     : ')
484            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_trp  : smv_i    : ')
485            DO jk = 1, nlay_i
486               CALL prt_ctl_info(' ')
487               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
488               CALL prt_ctl_info('   ~~~~~~~')
489               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp  : t_i      : ')
490               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp  : e_i      : ')
491            END DO
492         END DO
493      ENDIF
494      !
495      CALL wrk_dealloc( jpi,jpj,           zsm, zs0at, zatold, zeiold, zesold )
496      CALL wrk_dealloc( jpi,jpj,jpl,       zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e )
497      CALL wrk_dealloc( jpi,jpj,1,         zs0ow )
498      CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, zs0e )
499      CALL wrk_dealloc( jpi,jpj,jpl,       zviold, zvsold, zhimax )
500      !
501      IF( nn_timing == 1 )  CALL timing_stop('limtrp')
502   END SUBROUTINE lim_trp
503
504#else
505   !!----------------------------------------------------------------------
506   !!   Default option         Empty Module                No sea-ice model
507   !!----------------------------------------------------------------------
508CONTAINS
509   SUBROUTINE lim_trp        ! Empty routine
510   END SUBROUTINE lim_trp
511#endif
512
513   !!======================================================================
514END MODULE limtrp
Note: See TracBrowser for help on using the repository browser.