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

Last change on this file since 989 was 989, checked in by rblod, 16 years ago

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