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

source: branches/dev_002_LIM/NEMO/LIM_SRC_3/limtrp.F90 @ 830

Last change on this file since 830 was 825, checked in by ctlod, 16 years ago

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