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

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

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 3750

Last change on this file since 3750 was 3625, checked in by acc, 11 years ago

Branch dev_NOC_2012_r3555. #1006. Step 7. Check in code now merged with dev_r3385_NOCS04_HAMF

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