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 @ 888

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

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