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

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

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