New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limtrp.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta
  • Property svn:keywords set to Id
File size: 28.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'                                      LIM3 sea-ice model
9   !!----------------------------------------------------------------------
10   !!   lim_trp      : advection/diffusion process of sea ice
11   !!   lim_trp_init : initialization and namelist read
12   !!----------------------------------------------------------------------
13   USE phycst
14   USE dom_oce
15   USE in_out_manager  ! I/O manager
16   USE sbc_oce         ! Surface boundary condition: ocean fields
17   USE dom_ice
18   USE ice
19   USE limadv
20   USE limhdf
21   USE lbclnk
22   USE lib_mpp
23   USE par_ice
24   USE prtctl          ! Print control
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Routine accessibility
30   PUBLIC lim_trp       ! called by ice_step
31
32   !! * Shared module variables
33   REAL(wp), PUBLIC  ::   &  !:
34      bound  = 0.e0           !: boundary condit. (0.0 no-slip, 1.0 free-slip)
35
36   !! * Module variables
37   REAL(wp)  ::           &  ! constant values
38      epsi06 = 1.e-06  ,  &
39      epsi03 = 1.e-03  ,  &
40      epsi16 = 1.e-16  ,  &
41      rzero  = 0.e0    ,  &
42      rone   = 1.e0    ,  &
43      zeps10 = 1.e-10
44
45   !! * Substitution
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
49   !! $Id$
50   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52
53CONTAINS
54
55   SUBROUTINE lim_trp( kt ) 
56      !!-------------------------------------------------------------------
57      !!                   ***  ROUTINE lim_trp ***
58      !!                   
59      !! ** purpose : advection/diffusion process of sea ice
60      !!
61      !! ** method  : variables included in the process are scalar,   
62      !!     other values are considered as second order.
63      !!     For advection, a second order Prather scheme is used. 
64      !!
65      !! ** action :
66      !!
67      !! History :
68      !!   1.0  !  00-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet)  Original code
69      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
70      !!   2.0  !  04-01 (G. Madec, C. Ethe)  F90, mpp
71      !!   3.0  !  05-11 (M. Vancoppenolle)   Multi-layer sea ice, salinity variations
72      !!---------------------------------------------------------------------
73      INTEGER, INTENT(in) ::   kt     ! number of iteration
74      !! * Local Variables
75      INTEGER  ::   ji, jj, jk, jl, layer, &  ! dummy loop indices
76         initad           ! number of sub-timestep for the advection
77      INTEGER  ::   ji_maxu, ji_maxv, jj_maxu, jj_maxv
78
79      REAL(wp) ::  &                             
80         zindb  ,  &
81         zindsn ,  &
82         zindic ,  &
83         zusvosn,  &
84         zusvoic,  &
85         zvbord ,  &
86         zcfl   ,  &
87         zusnit ,  &
88         zrtt, zsal, zage, &
89         zbigval, ze, &
90         zmaxu, zmaxv
91
92      REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace
93         zui_u , zvi_v , zsm   ,         &
94         zs0at, zs0ow
95
96      REAL(wp), DIMENSION(jpi,jpj,jpl):: &  ! temporary workspace
97         zs0ice, zs0sn, zs0a   ,         &
98         zs0c0 ,                         &
99         zs0sm , zs0oi
100
101      ! MHE Multilayer heat content
102      REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl)  ::   &  ! temporary workspace
103         zs0e
104
105      !---------------------------------------------------------------------
106
107      IF( numit == nstart  )   CALL lim_trp_init      ! Initialization (first time-step only)
108
109      zsm(:,:) = area(:,:)
110
111      IF( ln_limdyn ) THEN
112         IF( kt == nit000 .AND. lwp ) THEN
113            WRITE(numout,*) ' lim_trp : Ice Advection'
114            WRITE(numout,*) ' ~~~~~~~'
115         ENDIF
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, fs_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( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN        !==  odd ice time step:  adv_x then adv_y  ==!
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(:,:), sxyopw(:,:) )
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_m(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      IF(ln_ctl) THEN   ! Control print
555         CALL prt_ctl_info(' ')
556         CALL prt_ctl_info(' - Cell values : ')
557         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
558         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp  : cell area :')
559         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp  : at_i      :')
560         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp  : vt_i      :')
561         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp  : vt_s      :')
562         DO jl = 1, jpl
563            CALL prt_ctl_info(' ')
564            CALL prt_ctl_info(' - Category : ', ivar1=jl)
565            CALL prt_ctl_info('   ~~~~~~~~~~')
566            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_trp  : a_i      : ')
567            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_trp  : ht_i     : ')
568            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_trp  : ht_s     : ')
569            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_trp  : v_i      : ')
570            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_trp  : v_s      : ')
571            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_trp  : e_s      : ')
572            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_trp  : t_su     : ')
573            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_trp  : t_snow   : ')
574            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_trp  : sm_i     : ')
575            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_trp  : smv_i    : ')
576            DO jk = 1, nlay_i
577               CALL prt_ctl_info(' ')
578               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
579               CALL prt_ctl_info('   ~~~~~~~')
580               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp  : t_i      : ')
581               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp  : e_i      : ')
582            END DO
583         END DO
584      ENDIF
585
586   END SUBROUTINE lim_trp
587
588
589   SUBROUTINE lim_trp_init
590      !!-------------------------------------------------------------------
591      !!                  ***  ROUTINE lim_trp_init  ***
592      !!
593      !! ** Purpose :   initialization of ice advection parameters
594      !!
595      !! ** Method  : Read the namicetrp namelist and check the parameter
596      !!       values called at the first timestep (nit000)
597      !!
598      !! ** input   :   Namelist namicetrp
599      !!
600      !! history :
601      !!   2.0  !  03-08 (C. Ethe)  Original code
602      !!-------------------------------------------------------------------
603      NAMELIST/namicetrp/ bound
604      !!-------------------------------------------------------------------
605
606      ! Read Namelist namicetrp
607      REWIND ( numnam_ice )
608      READ   ( numnam_ice  , namicetrp )
609      IF(lwp) THEN
610         WRITE(numout,*)
611         WRITE(numout,*) 'lim_trp_init : Ice parameters for advection '
612         WRITE(numout,*) '~~~~~~~~~~~~'
613         WRITE(numout,*) ' boundary conditions (0.0 no-slip, 1.0 free-slip) bound  = ', bound
614         WRITE(numout,*) 
615      ENDIF
616
617   END SUBROUTINE lim_trp_init
618
619#else
620   !!----------------------------------------------------------------------
621   !!   Default option         Empty Module                No sea-ice model
622   !!----------------------------------------------------------------------
623CONTAINS
624   SUBROUTINE lim_trp        ! Empty routine
625   END SUBROUTINE lim_trp
626#endif
627
628   !!======================================================================
629END MODULE limtrp
Note: See TracBrowser for help on using the repository browser.