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 trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 5134

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

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