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

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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 13 years ago

update licence of all NEMO files...

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