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

source: trunk/NEMO/LIM_SRC_3/limtrp.F90 @ 862

Last change on this file since 862 was 834, checked in by ctlod, 16 years ago

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