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/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 2833

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

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