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 on Ticket #154 – Attachment – NEMO

Ticket #154: limtrp.F90

File limtrp.F90, 28.1 KB (added by nemo_user, 16 years ago)
Line 
1MODULE limtrp
2   !!======================================================================
3   !!                       ***  MODULE limtrp   ***
4   !! LIM transport ice model : sea-ice advection/diffusion
5   !!======================================================================
6#if defined key_lim3
7   !!----------------------------------------------------------------------
8   !!   'key_lim3'                                      LIM3 sea-ice model
9   !!----------------------------------------------------------------------
10   !!   lim_trp      : advection/diffusion process of sea ice
11   !!   lim_trp_init : initialization and namelist read
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE phycst
15   USE dom_oce
16   USE daymod
17   USE in_out_manager  ! I/O manager
18   USE ice_oce         ! ice variables
19   USE sbc_oce         ! Surface boundary condition: ocean fields
20   USE dom_ice
21   USE ice
22   USE iceini
23   USE limistate
24   USE limadv
25   USE limhdf
26   USE lbclnk
27   USE lib_mpp
28   USE par_ice
29   USE prtctl          ! Print control
30
31   IMPLICIT NONE
32   PRIVATE
33
34   !! * Routine accessibility
35   PUBLIC lim_trp       ! called by ice_step
36
37   !! * Shared module variables
38   REAL(wp), PUBLIC  ::   &  !:
39      bound  = 0.e0           !: boundary condit. (0.0 no-slip, 1.0 free-slip)
40
41   !! * Module variables
42   REAL(wp)  ::           &  ! constant values
43      epsi06 = 1.e-06  ,  &
44      epsi03 = 1.e-03  ,  &
45      epsi16 = 1.e-16  ,  &
46      rzero  = 0.e0    ,  &
47      rone   = 1.e0    ,  &
48      zeps10 = 1.e-10
49
50   !! * Substitution
51#  include "vectopt_loop_substitute.h90"
52   !!----------------------------------------------------------------------
53   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
54   !! $ Id: $
55   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE lim_trp( kt ) 
60      !!-------------------------------------------------------------------
61      !!                   ***  ROUTINE lim_trp ***
62      !!                   
63      !! ** purpose : advection/diffusion process of sea ice
64      !!
65      !! ** method  : variables included in the process are scalar,   
66      !!     other values are considered as second order.
67      !!     For advection, a second order Prather scheme is used. 
68      !!
69      !! ** action :
70      !!
71      !! History :
72      !!   1.0  !  00-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet)  Original code
73      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
74      !!   2.0  !  04-01 (G. Madec, C. Ethe)  F90, mpp
75      !!   3.0  !  05-11 (M. Vancoppenolle)   Multi-layer sea ice, salinity variations
76      !!---------------------------------------------------------------------
77      INTEGER, INTENT(in) ::   kt     ! number of iteration
78      !! * Local Variables
79      INTEGER  ::   ji, jj, jk, jl, layer, &  ! dummy loop indices
80         initad           ! number of sub-timestep for the advection
81      INTEGER  ::   ji_maxu, ji_maxv, jj_maxu, jj_maxv
82      INTEGER  ::   kk
83
84      REAL(wp) ::  &                             
85         zindb  ,  &
86         zindsn ,  &
87         zindic ,  &
88         zusvosn,  &
89         zusvoic,  &
90         zvbord ,  &
91         zcfl   ,  &
92         zusnit ,  &
93         zrtt, zsal, zage, &
94         zbigval, ze, &
95         zmaxu, zmaxv
96
97      REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace
98         zui_u , zvi_v , zsm   ,         &
99         zs0at, zs0ow
100
101      REAL(wp), DIMENSION(jpi,jpj,jpl):: &  ! temporary workspace
102         zs0ice, zs0sn, zs0a   ,         &
103         zs0c0 ,                         &
104         zs0sm , zs0oi
105
106      ! MHE Multilayer heat content
107      REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl)  ::   &  ! temporary workspace
108         zs0e
109
110      !---------------------------------------------------------------------
111
112      IF( numit == nstart  )   CALL lim_trp_init      ! Initialization (first time-step only)
113
114      zsm(:,:) = area(:,:)
115
116      IF( ln_limdyn ) THEN
117         IF( kt == nit000  .AND. lwp) THEN
118            WRITE(numout,*) ' lim_trp : Ice Advection'
119            WRITE(numout,*) ' ~~~~~~~'
120         ENDIF
121
122         !-----------------------------------------------------------------------------!
123         ! 1) CFL Test                                                             
124         !-----------------------------------------------------------------------------!
125
126         !------------------------------------------
127         ! ice velocities at ocean U- and V-points
128         !------------------------------------------
129
130         ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions.       
131         zvbord = 1.0 + ( 1.0 - bound )
132         DO jj = 1, jpjm1
133            DO ji = 1, fs_jpim1
134               zui_u(ji,jj) = u_ice(ji,jj)
135               zvi_v(ji,jj) = v_ice(ji,jj)
136            END DO
137         END DO
138
139         ! Lateral boundary conditions
140         CALL lbc_lnk( zui_u, 'U', -1. )
141         CALL lbc_lnk( zvi_v, 'V', -1. )
142
143         !-------------------------
144         ! CFL test for stability
145         !-------------------------
146
147         zcfl  = 0.e0
148         zcfl  = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, :     ) ) * rdt_ice / e1u(1:jpim1, :     ) ) )
149         zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) )
150
151         zmaxu = 0.0
152         zmaxv = 0.0
153         DO ji = 1, jpim1
154            DO jj = 1, jpjm1
155               IF ( (ABS(zui_u(ji,jj)) .GT. zmaxu) ) THEN
156                  zmaxu = MAX(zui_u(ji,jj), zmaxu )
157                  ji_maxu = ji
158                  jj_maxu = jj
159               ENDIF
160               IF ( (ABS(zvi_v(ji,jj)) .GT. zmaxv) ) THEN
161                  zmaxv = MAX(zvi_v(ji,jj), zmaxv )
162                  ji_maxv = ji
163                  jj_maxv = jj
164               ENDIF
165            END DO
166         END DO
167
168         IF (lk_mpp ) CALL mpp_max(zcfl)
169
170         IF ( zcfl > 0.5 .AND. lwp ) &
171            WRITE(numout,*) 'lim_trp : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl
172
173         !-----------------------------------------------------------------------------!
174         ! 2) Computation of transported fields                                       
175         !-----------------------------------------------------------------------------!
176
177         !------------------------------------------------------
178         ! 1.1) Snow vol, ice vol, salt and age contents, area
179         !------------------------------------------------------
180
181         zs0ow (:,:) =  ato_i(:,:)    * area(:,:)           ! Open water area
182         DO jl = 1, jpl  !sum over thickness categories
183            ! area -> is the unmasked and masked area of T-grid cell
184            zs0sn (:,:,jl) =  v_s(:,:,jl)    * area(:,:)    ! Snow volume.
185            zs0ice(:,:,jl) =  v_i(:,:,jl)    * area(:,:)    ! Ice volume.
186            zs0a  (:,:,jl) =  a_i(:,:,jl)    * area(:,:)    ! Ice area
187            zs0sm (:,:,jl) =  smv_i(:,:,jl)  * area(:,:)    ! Salt content
188            zs0oi (:,:,jl) =  oa_i (:,:,jl)  * area(:,:)    ! Age content
189
190            !----------------------------------
191            ! 1.2) Ice and snow heat contents
192            !----------------------------------
193
194            zs0c0 (:,:,jl)     = e_s(:,:,1,jl)              ! Snow heat cont.
195            DO jk = 1, nlay_i
196               zs0e(:,:,jk,jl) = e_i(:,:,jk,jl)             ! Ice heat content
197            END DO
198         END DO
199
200         !-----------------------------------------------------------------------------!
201         ! 3) Advection of Ice fields                                             
202         !-----------------------------------------------------------------------------!
203
204         ! If ice drift field is too fast, use an appropriate time step for advection.         
205         initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )
206         zusnit = 1.0 / REAL( initad ) 
207
208         IF ( MOD( nday , 2 ) == 0) THEN
209            DO jk = 1,initad
210               !--- ice open water area
211               CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0ow(:,:) , sxopw(:,:) , & 
212                  sxxopw(:,:), syopw(:,:) , & 
213                  syyopw(:,:), sxopw(:,:) )
214               CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0ow(:,:) , sxopw (:,:) , &
215                  sxxopw(:,:), syopw (:,:) , & 
216                  syyopw(:,:), sxyopw(:,:) )
217               DO jl = 1, jpl
218                  !--- ice volume  ---
219                  CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , & 
220                     sxxice(:,:,jl) , syice (:,:,jl) , & 
221                     syyice(:,:,jl) , sxyice(:,:,jl) )
222                  CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , &
223                     sxxice(:,:,jl) , syice (:,:,jl) , & 
224                     syyice(:,:,jl) , sxyice(:,:,jl) )
225                  !--- snow volume  ---
226                  CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , &
227                     sxxsn (:,:,jl) , sysn  (:,:,jl) , &
228                     syysn (:,:,jl) , sxysn (:,:,jl) )
229                  CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , &
230                     sxxsn (:,:,jl) , sysn  (:,:,jl) , &
231                     syysn (:,:,jl) , sxysn (:,:,jl) )
232                  !--- ice salinity ---
233                  CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
234                     sxxsal(:,:,jl) , sysal (:,:,jl) , &
235                     syysal(:,:,jl) , sxysal(:,:,jl)  )
236                  CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
237                     sxxsal(:,:,jl) , sysal (:,:,jl) , &
238                     syysal(:,:,jl) , sxysal(:,:,jl)  )
239                  !--- ice age      ---     
240                  CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
241                     sxxage(:,:,jl) , syage (:,:,jl) , &
242                     syyage(:,:,jl) , sxyage(:,:,jl)  )
243                  CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
244                     sxxage(:,:,jl) , syage (:,:,jl) , &
245                     syyage(:,:,jl) , sxyage(:,:,jl)  )
246                  !--- ice concentrations ---
247                  CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , &
248                     sxxa  (:,:,jl) , sya   (:,:,jl) , & 
249                     syya  (:,:,jl) , sxya  (:,:,jl)  )
250                  CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , & 
251                     sxxa  (:,:,jl) , sya   (:,:,jl) , & 
252                     syya  (:,:,jl) , sxya  (:,:,jl)  )
253                  !--- ice / snow thermal energetic contents ---
254                  CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , & 
255                     sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
256                     syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
257                  CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , &
258                     sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
259                     syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
260                  DO layer = 1, nlay_i
261                     CALL lim_adv_x( zusnit, zui_u, rone , zsm, &
262                        zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , & 
263                        sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , &
264                        syye(:,:,layer,jl) , sxye(:,:,layer,jl) )
265                     CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, & 
266                        zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , & 
267                        sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , &
268                        syye(:,:,layer,jl) , sxye(:,:,layer,jl) )
269                  END DO
270               END DO
271            END DO
272         ELSE
273            DO jk = 1, initad
274               !--- ice volume  ---
275               CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0ow (:,:) , sxopw (:,:) , &
276                  sxxopw(:,:) , syopw (:,:) , & 
277                  syyopw(:,:) , sxyopw(:,:) )
278               CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0ow (:,:) , sxopw (:,:) , & 
279                  sxxopw(:,:) , syopw (:,:) , &
280                  syyopw(:,:) , sxyopw(:,:) )
281               DO jl = 1, jpl
282                  !--- ice volume  ---
283                  CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , &
284                     sxxice(:,:,jl) , syice (:,:,jl) , & 
285                     syyice(:,:,jl) , sxyice(:,:,jl) )
286                  CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , & 
287                     sxxice(:,:,jl) , syice (:,:,jl) , &
288                     syyice(:,:,jl) , sxyice(:,:,jl) )
289                  !--- snow volume  ---
290                  CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , & 
291                     sxxsn (:,:,jl) , sysn  (:,:,jl) , & 
292                     syysn (:,:,jl) , sxysn (:,:,jl) )
293                  CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , & 
294                     sxxsn (:,:,jl) , sysn  (:,:,jl) , & 
295                     syysn (:,:,jl) , sxysn (:,:,jl) )
296                  !--- ice salinity ---
297                  CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
298                     sxxsal(:,:,jl) , sysal (:,:,jl) , &
299                     syysal(:,:,jl) , sxysal(:,:,jl) )
300                  CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
301                     sxxsal(:,:,jl) , sysal (:,:,jl) , &
302                     syysal(:,:,jl) , sxysal(:,:,jl) )
303                  !--- ice age      ---
304                  CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
305                     sxxage(:,:,jl) , syage (:,:,jl) , & 
306                     syyage(:,:,jl) , sxyage(:,:,jl)  )
307                  CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
308                     sxxage(:,:,jl) , syage (:,:,jl) , &
309                     syyage(:,:,jl) , sxyage(:,:,jl)   )
310                  !--- ice concentration ---
311                  CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , & 
312                     sxxa  (:,:,jl) , sya   (:,:,jl) , & 
313                     syya  (:,:,jl) , sxya  (:,:,jl)  )
314                  CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , & 
315                     sxxa  (:,:,jl) , sya   (:,:,jl) , & 
316                     syya  (:,:,jl) , sxya  (:,:,jl)  )
317                  !--- ice / snow thermal energetic contents ---
318                  CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , & 
319                     sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
320                     syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
321                  CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , &
322                     sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
323                     syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
324                  DO layer = 1, nlay_i
325                     CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0e(:,:,layer,jl) , &
326                        sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , &
327                        syye (:,:,layer,jl), sxye (:,:,layer,jl) )
328                     CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0e(:,:,layer,jl) , &
329                        sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , &
330                        syye (:,:,layer,jl), sxye (:,:,layer,jl)  )
331                  END DO
332
333               END DO
334            END DO
335         ENDIF
336
337         !-------------------------------------------
338         ! Recover the properties from their contents
339         !-------------------------------------------
340
341         zs0ow (:,:)       = zs0ow(:,:) / area(:,:)
342         DO jl = 1, jpl
343            zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:)
344            zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:)
345            zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:)
346            zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:)
347            zs0a  (:,:,jl) = zs0a  (:,:,jl) / area(:,:)
348            zs0c0 (:,:,jl) = zs0c0 (:,:,jl) / area(:,:)
349            DO jk = 1, nlay_i
350               zs0e(:,:,jk,jl) = zs0e(:,:,jk,jl) / area(:,:)
351            END DO
352         END DO
353
354         !------------------------------------------------------------------------------!
355         ! 4) Diffusion of Ice fields                 
356         !------------------------------------------------------------------------------!
357
358         !------------------------------------
359         ! 4.1) diffusion of open water area
360         !------------------------------------
361
362         ! Compute total ice fraction
363         zs0at(:,:) = 0.0
364         DO jl = 1, jpl
365            DO jj = 1, jpj
366               DO ji = 1, jpi
367                  zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) !
368               END DO
369            END DO
370         END DO
371
372         ! Masked eddy diffusivity coefficient at ocean U- and V-points
373         DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
374            DO ji = 1 , fs_jpim1   ! vector opt.
375               pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji  ,jj) ) ) )   &
376                  &        * ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)
377               pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0at(ji,jj  ) ) ) )   &
378                  &        * ( 1.0 - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj)
379            END DO !jj
380         END DO ! ji
381
382         ! Diffusion
383         CALL lim_hdf( zs0ow (:,:) )
384
385         !----------------------------------------
386         ! 4.2) Diffusion of other ice variables
387         !----------------------------------------
388         DO jl = 1, jpl
389
390            ! Masked eddy diffusivity coefficient at ocean U- and V-points
391            DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
392               DO ji = 1 , fs_jpim1   ! vector opt.
393                  pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji  ,jj,jl) ) ) )   &
394                     &        * ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
395                  pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji,jj,jl  ) ) ) )   &
396                     &        * ( 1.0 - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
397               END DO !jj
398            END DO ! ji
399
400            CALL lim_hdf( zs0ice (:,:,jl) )
401            CALL lim_hdf( zs0sn  (:,:,jl) )
402            CALL lim_hdf( zs0sm  (:,:,jl) )
403            CALL lim_hdf( zs0oi  (:,:,jl) )
404            CALL lim_hdf( zs0a   (:,:,jl) )
405            CALL lim_hdf( zs0c0  (:,:,jl) )
406            DO jk = 1, nlay_i
407               CALL lim_hdf( zs0e (:,:,jk,jl) )
408            END DO ! jk
409         END DO !jl
410
411         !-----------------------------------------
412         ! 4.3) Remultiply everything by ice area
413         !-----------------------------------------
414         zs0ow(:,:) = MAX(rzero, zs0ow(:,:) * area(:,:) )
415         DO jl = 1, jpl
416            zs0ice(:,:,jl) = MAX( rzero, zs0ice(:,:,jl) * area(:,:) )    !!bug:  est-ce utile
417            zs0sn (:,:,jl) = MAX( rzero, zs0sn (:,:,jl) * area(:,:) )    !!bug:  cf /area  juste apres
418            zs0sm (:,:,jl) = MAX( rzero, zs0sm (:,:,jl) * area(:,:) )    !!bug:  cf /area  juste apres
419            zs0oi (:,:,jl) = MAX( rzero, zs0oi (:,:,jl) * area(:,:) )
420            zs0a  (:,:,jl) = MAX( rzero, zs0a  (:,:,jl) * area(:,:) )    !! suppress both change le resultat
421            zs0c0 (:,:,jl) = MAX( rzero, zs0c0 (:,:,jl) * area(:,:) )
422            DO jk = 1, nlay_i
423               zs0e(:,:,jk,jl) = MAX( rzero, zs0e (:,:,jk,jl) * area(:,:) )
424            END DO ! jk
425         END DO ! jl
426
427         !------------------------------------------------------------------------------!
428         ! 5) Update and limit ice properties after transport                           
429         !------------------------------------------------------------------------------!
430
431         !--------------------------------------------------
432         ! 5.1) Recover mean values over the grid squares.
433         !--------------------------------------------------
434
435         DO jl = 1, jpl
436            DO jk = 1, nlay_i
437               DO jj = 1, jpj
438                  DO ji = 1, jpi
439                     zs0e (ji,jj,jk,jl) =  &
440                        MAX( rzero, zs0e (ji,jj,jk,jl) / area(ji,jj) )
441                  END DO
442               END DO
443            END DO
444         END DO
445
446         DO jj = 1, jpj
447            DO ji = 1, jpi
448               zs0ow (ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) )
449            END DO
450         END DO
451
452         zs0at(:,:) = 0.0
453         DO jl = 1, jpl
454            DO jj = 1, jpj
455               DO ji = 1, jpi
456                  zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl)/area(ji,jj) )
457                  zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl)/area(ji,jj) )
458                  zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl)/area(ji,jj) )
459                  zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl)/area(ji,jj) )
460                  zs0a  (ji,jj,jl) = MAX( rzero, zs0a  (ji,jj,jl)/area(ji,jj) )
461                  zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl)/area(ji,jj) )
462                  zs0at (ji,jj)    = zs0at(ji,jj) + zs0a(ji,jj,jl)
463               END DO
464            END DO
465         END DO
466
467         !---------------------------------------------------------
468         ! 5.2) Snow thickness, Ice thickness, Ice concentrations
469         !---------------------------------------------------------
470
471         DO jj = 1, jpj
472            DO ji = 1, jpi
473               zindb         = MAX( 0.0 , SIGN( 1.0, zs0at(ji,jj) - zeps10) )
474               zs0ow(ji,jj)  = (1.0 - zindb) + zindb*MAX( zs0ow(ji,jj), 0.00 )
475               ato_i(ji,jj)  = zs0ow(ji,jj)
476            END DO
477         END DO
478
479         ! Remove very small areas
480         DO jl = 1, jpl
481            DO jj = 1, jpj
482               DO ji = 1, jpi
483                  zindb         = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - zeps10) )
484
485                  zs0a(ji,jj,jl)  = zindb * MIN( zs0a(ji,jj,jl), 0.99 )
486                  v_s(ji,jj,jl)   = zindb * zs0sn (ji,jj,jl) 
487                  v_i(ji,jj,jl)   = zindb * zs0ice(ji,jj,jl)
488
489                  zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )
490                  zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
491                  zindb           = MAX( zindsn, zindic )
492                  zs0a (ji,jj,jl) = zindb  * zs0a(ji,jj,jl) !ice concentration
493                  a_i  (ji,jj,jl) = zs0a(ji,jj,jl)
494                  v_s  (ji,jj,jl) = zindsn * v_s(ji,jj,jl)
495                  v_i  (ji,jj,jl) = zindic * v_i(ji,jj,jl)
496               END DO
497            END DO
498         END DO
499
500         DO jj = 1, jpj
501            DO ji = 1, jpi
502               zs0at(ji,jj)       = SUM( zs0a(ji,jj,1:jpl) )
503            END DO
504         END DO
505
506         !----------------------
507         ! 5.3) Ice properties
508         !----------------------
509
510         zbigval         =  1.0d+13
511
512         DO jl = 1, jpl
513            DO jj = 1, jpj
514               DO ji = 1, jpi
515
516                  ! Switches and dummy variables
517                  zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi16 )
518                  zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi16 )
519                  zrtt            = 173.15 * rone 
520                  zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )
521                  zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
522                  zindb           = MAX( zindsn, zindic )
523
524                  ! Ice salinity and age
525                  zsal            = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj)  , &
526                     zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * &
527                     v_i(ji,jj,jl)
528                  IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 
529                     smv_i(ji,jj,jl) = zindic*zsal + (1.0-zindic)*0.0
530
531                  zage            = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / & 
532                     MAX( a_i(ji,jj,jl), epsi16 )  ), 0.0 ) * &
533                     a_i(ji,jj,jl)
534                  oa_i (ji,jj,jl)  = zindic*zage 
535
536                  ! Snow heat content
537                  ze              =  MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval )
538                  e_s(ji,jj,1,jl) = zindsn * ze + (1.0 - zindsn) * 0.0     
539
540               END DO !ji
541            END DO !jj
542         END DO ! jl
543
544         DO jl = 1, jpl
545            DO jk = 1, nlay_i
546               DO jj = 1, jpj
547                  DO ji = 1, jpi
548                     ! Ice heat content
549                     zindic          =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
550                     ze              =  MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval )
551                     e_i(ji,jj,jk,jl) = zindic * ze    + ( 1.0 - zindic ) * 0.0
552                  END DO !ji
553               END DO ! jj
554            END DO ! jk
555         END DO ! jl
556
557      ENDIF
558
559      IF(ln_ctl) THEN   ! Control print
560         CALL prt_ctl_info(' ')
561         CALL prt_ctl_info(' - Cell values : ')
562         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
563         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp  : cell area :')
564         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp  : at_i      :')
565         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp  : vt_i      :')
566         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp  : vt_s      :')
567         DO jl = 1, jpl
568            CALL prt_ctl_info(' ')
569            CALL prt_ctl_info(' - Category : ', ivar1=jl)
570            CALL prt_ctl_info('   ~~~~~~~~~~')
571            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_trp  : a_i      : ')
572            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_trp  : ht_i     : ')
573            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_trp  : ht_s     : ')
574            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_trp  : v_i      : ')
575            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_trp  : v_s      : ')
576            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_trp  : e_s      : ')
577            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_trp  : t_su     : ')
578            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_trp  : t_snow   : ')
579            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_trp  : sm_i     : ')
580            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_trp  : smv_i    : ')
581            DO jk = 1, nlay_i
582               CALL prt_ctl_info(' ')
583               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
584               CALL prt_ctl_info('   ~~~~~~~')
585               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp  : t_i      : ')
586               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp  : e_i      : ')
587            END DO
588         END DO
589      ENDIF
590
591   END SUBROUTINE lim_trp
592
593
594   SUBROUTINE lim_trp_init
595      !!-------------------------------------------------------------------
596      !!                  ***  ROUTINE lim_trp_init  ***
597      !!
598      !! ** Purpose :   initialization of ice advection parameters
599      !!
600      !! ** Method  : Read the namicetrp namelist and check the parameter
601      !!       values called at the first timestep (nit000)
602      !!
603      !! ** input   :   Namelist namicetrp
604      !!
605      !! history :
606      !!   2.0  !  03-08 (C. Ethe)  Original code
607      !!-------------------------------------------------------------------
608      NAMELIST/namicetrp/ bound
609      !!-------------------------------------------------------------------
610
611      ! Read Namelist namicetrp
612      REWIND ( numnam_ice )
613      READ   ( numnam_ice  , namicetrp )
614      IF(lwp) THEN
615         WRITE(numout,*)
616         WRITE(numout,*) 'lim_trp_init : Ice parameters for advection '
617         WRITE(numout,*) '~~~~~~~~~~~~'
618         WRITE(numout,*) ' boundary conditions (0.0 no-slip, 1.0 free-slip) bound  = ', bound
619         WRITE(numout,*) 
620      ENDIF
621
622   END SUBROUTINE lim_trp_init
623
624#else
625   !!----------------------------------------------------------------------
626   !!   Default option         Empty Module                No sea-ice model
627   !!----------------------------------------------------------------------
628CONTAINS
629   SUBROUTINE lim_trp        ! Empty routine
630   END SUBROUTINE lim_trp
631#endif
632
633   !!======================================================================
634END MODULE limtrp