source: branches/2012/dev_INGV/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 3577

Last change on this file since 3577 was 3577, checked in by adani, 8 years ago

ticket #998. Step 5: Add in changes from the trunk between revisions 3521 and 3555.

  • Property svn:keywords set to Id
File size: 26.1 KB
Line 
1MODULE limtrp
2   !!======================================================================
3   !!                       ***  MODULE limtrp   ***
4   !! LIM transport ice model : sea-ice advection/diffusion
5   !!======================================================================
6   !! History : LIM-2 ! 2000-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet)  Original code
7   !!            3.0  ! 2005-11 (M. Vancoppenolle)   Multi-layer sea ice, salinity variations
8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3'                                      LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_trp      : advection/diffusion process of sea ice
15   !!----------------------------------------------------------------------
16   USE phycst          ! physical constant
17   USE dom_oce         ! ocean domain
18   USE sbc_oce         ! ocean surface boundary condition
19   USE par_ice         ! LIM-3 parameter
20   USE dom_ice         ! LIM-3 domain
21   USE ice             ! LIM-3 variables
22   USE limadv          ! LIM-3 advection
23   USE limhdf          ! LIM-3 horizontal diffusion
24   USE in_out_manager  ! I/O manager
25   USE lbclnk          ! lateral boundary conditions -- MPP exchanges
26   USE lib_mpp         ! MPP library
27   USE wrk_nemo        ! work arrays
28   USE prtctl          ! Print control
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   lim_trp    ! called by ice_step
34
35   REAL(wp)  ::   epsi06 = 1.e-06_wp   ! constant values
36   REAL(wp)  ::   epsi03 = 1.e-03_wp 
37   REAL(wp)  ::   zeps10 = 1.e-10_wp 
38   REAL(wp)  ::   epsi16 = 1.e-16_wp
39   REAL(wp)  ::   rzero  = 0._wp   
40   REAL(wp)  ::   rone   = 1._wp
41
42   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) ::   zs0e
43
44   !! * Substitution
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
48   !! $Id$
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE lim_trp( kt ) 
54      !!-------------------------------------------------------------------
55      !!                   ***  ROUTINE lim_trp ***
56      !!                   
57      !! ** purpose : advection/diffusion process of sea ice
58      !!
59      !! ** method  : variables included in the process are scalar,   
60      !!     other values are considered as second order.
61      !!     For advection, a second order Prather scheme is used. 
62      !!
63      !! ** action :
64      !!---------------------------------------------------------------------
65      INTEGER, INTENT(in) ::   kt   ! number of iteration
66      !
67      INTEGER  ::   ji, jj, jk, jl, layer   ! dummy loop indices
68      INTEGER  ::   initad                  ! number of sub-timestep for the advection
69      INTEGER  ::   ierr                    ! error status
70      REAL(wp) ::   zindb  , zindsn , zindic      ! local scalar
71      REAL(wp) ::   zusvosn, zusvoic, zbigval     !   -      -
72      REAL(wp) ::   zcfl , zusnit , zrtt          !   -      -
73      REAL(wp) ::   ze   , zsal   , zage          !   -      -
74      !
75      REAL(wp), POINTER, DIMENSION(:,:)      ::   zui_u, zvi_v, zsm, zs0at, zs0ow
76      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi
77      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e
78      !!---------------------------------------------------------------------
79
80      CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )
81      CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
82      CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e )
83
84      IF( numit == nstart .AND. lwp ) THEN
85         WRITE(numout,*)
86         IF( ln_limdyn ) THEN   ;   WRITE(numout,*) 'lim_trp : Ice transport '
87         ELSE                   ;   WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
88         ENDIF
89         WRITE(numout,*) '~~~~~~~~~~~~'
90      ENDIF
91     
92      zsm(:,:) = area(:,:)
93
94      !                             !-------------------------------------!
95      IF( ln_limdyn ) THEN          !   Advection of sea ice properties   !
96         !                          !-------------------------------------!
97         !
98
99         !-------------------------
100         ! transported fields                                       
101         !-------------------------
102         ! Snow vol, ice vol, salt and age contents, area
103         zs0ow(:,:) = ato_i(:,:) * area(:,:)               ! Open water area
104         DO jl = 1, jpl
105            zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume
106            zs0ice(:,:,jl)   = v_i  (:,:,jl) * area(:,:)    ! Ice  volume
107            zs0a  (:,:,jl)   = a_i  (:,:,jl) * area(:,:)    ! Ice area
108            zs0sm (:,:,jl)   = smv_i(:,:,jl) * area(:,:)    ! Salt content
109            zs0oi (:,:,jl)   = oa_i (:,:,jl) * area(:,:)    ! Age content
110            zs0c0 (:,:,jl)   = e_s  (:,:,1,jl)              ! Snow heat content
111            zs0e  (:,:,:,jl) = e_i  (:,:,:,jl)              ! Ice  heat content
112         END DO
113
114         !--------------------------
115         ! Advection of Ice fields  (Prather scheme)                                           
116         !--------------------------
117         ! If ice drift field is too fast, use an appropriate time step for advection.         
118         ! CFL test for stability
119         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) )
120         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) )
121         IF(lk_mpp )   CALL mpp_max( zcfl )
122!!gm more readability:
123!         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
124!         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
125!         ENDIF
126!!gm end
127         initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )
128         zusnit = 1.0 / REAL( initad ) 
129         IF( zcfl > 0.5 .AND. lwp )   &
130            WRITE(numout,*) 'lim_trp_2 : CFL violation at day ', nday, ', cfl = ', zcfl,   &
131               &                        ': the ice time stepping is split in two'
132
133         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
134            DO jk = 1,initad
135               CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area
136                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
137               CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:),   &
138                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
139               DO jl = 1, jpl
140                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
141                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
142                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
143                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
144                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
145                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
146                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
147                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
148                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
149                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
150                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
151                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
152                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---     
153                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
154                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &
155                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
156                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
157                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
158                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
159                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
160                  CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
161                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
162                  CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
163                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
164                  DO layer = 1, nlay_i                                                           !--- ice heat contents ---
165                     CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
166                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
167                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
168                     CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
169                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
170                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
171                  END DO
172               END DO
173            END DO
174         ELSE
175            DO jk = 1, initad
176               CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area
177                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
178               CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:),   &
179                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
180               DO jl = 1, jpl
181                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
182                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
183                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
184                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
185                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
186                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
187                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
188                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
189                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
190                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
191                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
192                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
193
194                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
195                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
196                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &
197                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
198                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
199                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
200                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &
201                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
202                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
203                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
204                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
205                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
206                  DO layer = 1, nlay_i                                                           !--- ice heat contents ---
207                     CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
208                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
209                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
210                     CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
211                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
212                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
213                  END DO
214               END DO
215            END DO
216         ENDIF
217
218         !-------------------------------------------
219         ! Recover the properties from their contents
220         !-------------------------------------------
221         zs0ow(:,:) = zs0ow(:,:) / area(:,:)
222         DO jl = 1, jpl
223            zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:)
224            zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:)
225            zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:)
226            zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:)
227            zs0a  (:,:,jl) = zs0a  (:,:,jl) / area(:,:)
228            zs0c0 (:,:,jl) = zs0c0 (:,:,jl) / area(:,:)
229            DO jk = 1, nlay_i
230               zs0e(:,:,jk,jl) = zs0e(:,:,jk,jl) / area(:,:)
231            END DO
232         END DO
233
234         !------------------------------------------------------------------------------!
235         ! 4) Diffusion of Ice fields                 
236         !------------------------------------------------------------------------------!
237
238         !--------------------------------
239         !  diffusion of open water area
240         !--------------------------------
241         zs0at(:,:) = zs0a(:,:,1)      ! total ice fraction
242         DO jl = 2, jpl
243            zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl)
244         END DO
245         !
246         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
247         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
248            DO ji = 1 , fs_jpim1   ! vector opt.
249               pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji  ,jj) ) ) )   &
250                  &        * ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)
251               pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0at(ji,jj  ) ) ) )   &
252                  &        * ( 1._wp - MAX( rzero, SIGN( rone,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj)
253            END DO
254         END DO
255         !
256         CALL lim_hdf( zs0ow (:,:) )   ! Diffusion
257
258         !------------------------------------
259         !  Diffusion of other ice variables
260         !------------------------------------
261         DO jl = 1, jpl
262         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
263            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
264               DO ji = 1 , fs_jpim1   ! vector opt.
265                  pahu(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji  ,jj,jl) ) ) )   &
266                     &        * ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
267                  pahv(ji,jj) = ( 1._wp - MAX( rzero, SIGN( rone, -zs0a(ji,jj  ,jl) ) ) )   &
268                     &        * ( 1._wp - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
269               END DO
270            END DO
271
272            CALL lim_hdf( zs0ice (:,:,jl) )
273            CALL lim_hdf( zs0sn  (:,:,jl) )
274            CALL lim_hdf( zs0sm  (:,:,jl) )
275            CALL lim_hdf( zs0oi  (:,:,jl) )
276            CALL lim_hdf( zs0a   (:,:,jl) )
277            CALL lim_hdf( zs0c0  (:,:,jl) )
278            DO jk = 1, nlay_i
279               CALL lim_hdf( zs0e (:,:,jk,jl) )
280            END DO
281         END DO
282
283         !-----------------------------------------
284         !  Remultiply everything by ice area
285         !-----------------------------------------
286         zs0ow(:,:) = MAX( rzero, zs0ow(:,:) * area(:,:) )
287         DO jl = 1, jpl
288            zs0ice(:,:,jl) = MAX( rzero, zs0ice(:,:,jl) * area(:,:) )    !!bug:  est-ce utile
289            zs0sn (:,:,jl) = MAX( rzero, zs0sn (:,:,jl) * area(:,:) )    !!bug:  cf /area  juste apres
290            zs0sm (:,:,jl) = MAX( rzero, zs0sm (:,:,jl) * area(:,:) )    !!bug:  cf /area  juste apres
291            zs0oi (:,:,jl) = MAX( rzero, zs0oi (:,:,jl) * area(:,:) )
292            zs0a  (:,:,jl) = MAX( rzero, zs0a  (:,:,jl) * area(:,:) )    !! suppress both change le resultat
293            zs0c0 (:,:,jl) = MAX( rzero, zs0c0 (:,:,jl) * area(:,:) )
294            DO jk = 1, nlay_i
295               zs0e(:,:,jk,jl) = MAX( rzero, zs0e (:,:,jk,jl) * area(:,:) )
296            END DO ! jk
297         END DO ! jl
298
299         !------------------------------------------------------------------------------!
300         ! 5) Update and limit ice properties after transport                           
301         !------------------------------------------------------------------------------!
302
303         !--------------------------------------------------
304         ! 5.1) Recover mean values over the grid squares.
305         !--------------------------------------------------
306
307         DO jl = 1, jpl
308            DO jk = 1, nlay_i
309               DO jj = 1, jpj
310                  DO ji = 1, jpi
311                     zs0e(ji,jj,jk,jl) = MAX( rzero, zs0e(ji,jj,jk,jl) / area(ji,jj) )
312                  END DO
313               END DO
314            END DO
315         END DO
316
317         DO jj = 1, jpj
318            DO ji = 1, jpi
319               zs0ow(ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) )
320            END DO
321         END DO
322
323         zs0at(:,:) = 0._wp
324         DO jl = 1, jpl
325            DO jj = 1, jpj
326               DO ji = 1, jpi
327                  zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl)/area(ji,jj) )
328                  zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl)/area(ji,jj) )
329                  zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl)/area(ji,jj) )
330                  zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl)/area(ji,jj) )
331                  zs0a  (ji,jj,jl) = MAX( rzero, zs0a  (ji,jj,jl)/area(ji,jj) )
332                  zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl)/area(ji,jj) )
333                  zs0at (ji,jj)    = zs0at(ji,jj) + zs0a(ji,jj,jl)
334               END DO
335            END DO
336         END DO
337
338         !---------------------------------------------------------
339         ! 5.2) Snow thickness, Ice thickness, Ice concentrations
340         !---------------------------------------------------------
341         DO jj = 1, jpj
342            DO ji = 1, jpi
343               zindb        = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - zeps10) )
344               zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp )
345               ato_i(ji,jj) = zs0ow(ji,jj)
346            END DO
347         END DO
348
349         DO jl = 1, jpl         ! Remove very small areas
350            DO jj = 1, jpj
351               DO ji = 1, jpi
352                  zindb         = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - zeps10) )
353                  !
354                  zs0a(ji,jj,jl) = zindb * MIN( zs0a(ji,jj,jl), 0.99 )
355                  v_s(ji,jj,jl)  = zindb * zs0sn (ji,jj,jl) 
356                  v_i(ji,jj,jl)  = zindb * zs0ice(ji,jj,jl)
357                  !
358                  zindsn         = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )
359                  zindic         = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
360                  zindb          = MAX( zindsn, zindic )
361                  zs0a(ji,jj,jl) = zindb  * zs0a(ji,jj,jl) !ice concentration
362                  a_i (ji,jj,jl) = zs0a(ji,jj,jl)
363                  v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl)
364                  v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl)
365               END DO
366            END DO
367         END DO
368
369         DO jj = 1, jpj
370            DO ji = 1, jpi
371               zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) )
372            END DO
373         END DO
374
375         !----------------------
376         ! 5.3) Ice properties
377         !----------------------
378
379         zbigval = 1.d+13
380
381         DO jl = 1, jpl
382            DO jj = 1, jpj
383               DO ji = 1, jpi
384
385                  ! Switches and dummy variables
386                  zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi16 )
387                  zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi16 )
388                  zrtt            = 173.15 * rone 
389                  zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )
390                  zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
391                  zindb           = MAX( zindsn, zindic )
392
393                  ! Ice salinity and age
394                  zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj)  , &
395                     zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl)
396                  IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 
397                     smv_i(ji,jj,jl) = zindic*zsal + (1.0-zindic)*0.0
398
399                  zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / & 
400                     MAX( a_i(ji,jj,jl), epsi16 )  ), 0.0 ) * a_i(ji,jj,jl)
401                  oa_i (ji,jj,jl)  = zindic*zage 
402
403                  ! Snow heat content
404                  ze              =  MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval )
405                  e_s(ji,jj,1,jl) = zindsn * ze + (1.0 - zindsn) * 0.0     
406
407               END DO !ji
408            END DO !jj
409         END DO ! jl
410
411         DO jl = 1, jpl
412            DO jk = 1, nlay_i
413               DO jj = 1, jpj
414                  DO ji = 1, jpi
415                     ! Ice heat content
416                     zindic          =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )
417                     ze              =  MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval )
418                     e_i(ji,jj,jk,jl) = zindic * ze    + ( 1.0 - zindic ) * 0.0
419                  END DO !ji
420               END DO ! jj
421            END DO ! jk
422         END DO ! jl
423
424      ENDIF
425
426      IF(ln_ctl) THEN   ! Control print
427         CALL prt_ctl_info(' ')
428         CALL prt_ctl_info(' - Cell values : ')
429         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
430         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp  : cell area :')
431         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp  : at_i      :')
432         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp  : vt_i      :')
433         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp  : vt_s      :')
434         DO jl = 1, jpl
435            CALL prt_ctl_info(' ')
436            CALL prt_ctl_info(' - Category : ', ivar1=jl)
437            CALL prt_ctl_info('   ~~~~~~~~~~')
438            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_trp  : a_i      : ')
439            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_trp  : ht_i     : ')
440            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_trp  : ht_s     : ')
441            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_trp  : v_i      : ')
442            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_trp  : v_s      : ')
443            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_trp  : e_s      : ')
444            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_trp  : t_su     : ')
445            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_trp  : t_snow   : ')
446            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_trp  : sm_i     : ')
447            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_trp  : smv_i    : ')
448            DO jk = 1, nlay_i
449               CALL prt_ctl_info(' ')
450               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
451               CALL prt_ctl_info('   ~~~~~~~')
452               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp  : t_i      : ')
453               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp  : e_i      : ')
454            END DO
455         END DO
456      ENDIF
457      !
458      CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow )
459      CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
460      CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e )
461      !
462   END SUBROUTINE lim_trp
463
464#else
465   !!----------------------------------------------------------------------
466   !!   Default option         Empty Module                No sea-ice model
467   !!----------------------------------------------------------------------
468CONTAINS
469   SUBROUTINE lim_trp        ! Empty routine
470   END SUBROUTINE lim_trp
471#endif
472
473   !!======================================================================
474END MODULE limtrp
Note: See TracBrowser for help on using the repository browser.