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/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 8321

Last change on this file since 8321 was 8321, checked in by clem, 7 years ago

STEP3 (1): clean separation between sea-ice and ocean
  • Property svn:keywords set to Id
File size: 29.5 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 , ONLY : nn_fsbc
19   USE ice            ! ice variables
20   USE limvar         !
21   USE limadv_prather ! advection scheme (Prather)
22   USE limadv_umx     ! advection scheme (ultimate-macho)
23   !
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   USE timing         ! Timing
31   USE limcons        ! conservation tests
32   USE limctl         ! control prints
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   lim_trp    ! called by icestp
38
39   INTEGER  ::   ncfl                 ! number of ice time step with CFL>1/2 
40
41   !! * Substitution
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE lim_trp( kt ) 
51      !!-------------------------------------------------------------------
52      !!                   ***  ROUTINE lim_trp ***
53      !!                   
54      !! ** purpose : advection/diffusion process of sea ice
55      !!
56      !! ** method  : variables included in the process are scalar,   
57      !!     other values are considered as second order.
58      !!     For advection, one can choose between
59      !!     a) an Ultimate-Macho scheme (whose order is defined by nn_limadv_ord) => nn_limadv=0
60      !!     b) and a second order Prather scheme => nn_limadv=-1
61      !!
62      !! ** action :
63      !!---------------------------------------------------------------------
64      INTEGER, INTENT(in) ::   kt   ! number of iteration
65      !
66      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
67      INTEGER  ::   initad                  ! number of sub-timestep for the advection
68      REAL(wp) ::   zcfl , zusnit           !   -      -
69      CHARACTER(len=80) :: cltmp
70      !
71      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
72      REAL(wp) ::    zdv, zda
73      REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold, zsmvold 
74      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax, zviold, zvsold
75      ! --- ultimate macho only --- !
76      REAL(wp)                               ::   zdt
77      REAL(wp), POINTER, DIMENSION(:,:)      ::   zudy, zvdx, zcu_box, zcv_box
78      ! --- prather only --- !
79      REAL(wp), POINTER, DIMENSION(:,:)      ::   zarea
80      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw
81      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
82      ! MV MP 2016
83      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ap , z0vp
84      REAL(wp) ::   za_old
85      ! END MV MP 2016
86      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei
87      !!
88      !!---------------------------------------------------------------------
89      IF( nn_timing == 1 )  CALL timing_start('limtrp')
90
91      CALL wrk_alloc( jpi,jpj,                            zatold, zeiold, zesold, zsmvold )
92      CALL wrk_alloc( jpi,jpj,jpl,                        zhimax, zviold, zvsold )
93 
94      IF( kt == nit000 .AND. lwp ) THEN
95         WRITE(numout,*)''
96         WRITE(numout,*)'limtrp'
97         WRITE(numout,*)'~~~~~~'
98         ncfl = 0                ! nb of time step with CFL > 1/2
99      ENDIF
100     
101      CALL lim_var_agg( 1 ) ! integrated values + ato_i
102
103      !-------------------------------------!
104      !   Advection of sea ice properties   !
105      !-------------------------------------!
106
107      ! conservation test
108      IF( ln_limdiachk )   CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
109     
110      ! store old values for diag
111      zviold = v_i
112      zvsold = v_s
113      zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 )
114      zeiold (:,:) = et_i
115      zesold (:,:) = et_s 
116
117      !--- Thickness correction init. --- !
118      zatold(:,:) = at_i
119      DO jl = 1, jpl
120         DO jj = 1, jpj
121            DO ji = 1, jpi
122               rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
123               ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
124               ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
125            END DO
126         END DO
127      END DO
128      ! --- Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick --- !
129      zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
130      DO jl = 1, jpl
131         DO jj = 2, jpjm1
132            DO ji = 2, jpim1
133               zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) )
134            END DO
135         END DO
136         CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
137      END DO
138         
139      ! --- If ice drift field is too fast, use an appropriate time step for advection --- !       
140      zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )    ! CFL test for stability
141      zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
142      IF( lk_mpp )   CALL mpp_max( zcfl )
143     
144      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
145      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
146      ENDIF
147     
148!!      IF( zcfl > 0.5_wp .AND. lwp ) THEN
149!!         ncfl = ncfl + 1
150!!         IF( ncfl > 0 ) THEN   
151!!            WRITE(cltmp,'(i6.1)') ncfl
152!!            CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
153!!         ENDIF
154!!      ENDIF
155
156      SELECT CASE ( nn_limadv )
157         
158                       !=============================!
159      CASE ( 0 )       !==  Ultimate-MACHO scheme  ==!                   
160                       !=============================!
161     
162         CALL wrk_alloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box )
163     
164         IF( kt == nit000 .AND. lwp ) THEN
165            WRITE(numout,*)''
166            WRITE(numout,*)'lim_adv_umx : Ultimate-MACHO advection scheme'
167            WRITE(numout,*)'~~~~~~~~~~~'
168         ENDIF
169         !
170         zdt = rdt_ice / REAL(initad)
171         
172         ! transport
173         zudy(:,:) = u_ice(:,:) * e2u(:,:)
174         zvdx(:,:) = v_ice(:,:) * e1v(:,:)
175         
176         ! define velocity for advection: u*grad(H)
177         DO jj = 2, jpjm1
178            DO ji = fs_2, fs_jpim1
179               IF    ( u_ice(ji,jj) * u_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
180               ELSEIF( u_ice(ji,jj)                  >  0._wp ) THEN   ;   zcu_box(ji,jj) = u_ice(ji-1,jj)
181               ELSE                                                    ;   zcu_box(ji,jj) = u_ice(ji  ,jj)
182               ENDIF
183               
184               IF    ( v_ice(ji,jj) * v_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
185               ELSEIF( v_ice(ji,jj)                  >  0._wp ) THEN   ;   zcv_box(ji,jj) = v_ice(ji,jj-1)
186               ELSE                                                    ;   zcv_box(ji,jj) = v_ice(ji,jj  )
187               ENDIF
188            END DO
189         END DO
190         
191         ! advection
192         DO jt = 1, initad
193            CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, ato_i(:,:) )       ! Open water area
194            DO jl = 1, jpl
195               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_i(:,:,jl) )      ! Ice area
196               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_i(:,:,jl) )      ! Ice  volume
197               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, smv_i(:,:,jl) )    ! Salt content
198               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, oa_i (:,:,jl) )    ! Age content
199               DO jk = 1, nlay_i
200                  CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_i(:,:,jk,jl) )   ! Ice  heat content
201               END DO
202               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_s(:,:,jl) )      ! Snow volume
203               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_s(:,:,1,jl) )    ! Snow heat content
204               ! MV MP 2016
205               IF ( nn_pnd_scheme > 0 ) THEN
206                  CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_ip(:,:,jl) )  ! Melt pond fraction
207                  CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_ip(:,:,jl) )  ! Melt pond volume
208               ENDIF
209               ! END MV MP 2016
210            END DO
211         END DO
212         !
213         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
214         DO jl = 2, jpl
215            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
216         END DO
217         !
218         CALL wrk_dealloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box )
219         
220                       !=============================!
221      CASE ( -1 )      !==     Prather scheme      ==!                   
222                       !=============================!
223
224         CALL wrk_alloc( jpi,jpj,            zarea )
225         CALL wrk_alloc( jpi,jpj,1,          z0opw )
226         CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
227         CALL wrk_alloc( jpi,jpj,jpl,        z0ap , z0vp )
228         CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
229         
230         IF( kt == nit000 .AND. lwp ) THEN
231            WRITE(numout,*)''
232            WRITE(numout,*)'lim_adv_xy : Prather advection scheme'
233            WRITE(numout,*)'~~~~~~~~~~~'
234         ENDIF
235         
236         zarea(:,:) = e1e2t(:,:)
237         
238         !-------------------------
239         ! transported fields                                       
240         !-------------------------
241         z0opw(:,:,1) = ato_i(:,:) * e1e2t(:,:)             ! Open water area
242         DO jl = 1, jpl
243            z0snw (:,:,jl)  = v_s  (:,:,  jl) * e1e2t(:,:)  ! Snow volume
244            z0ice(:,:,jl)   = v_i  (:,:,  jl) * e1e2t(:,:)  ! Ice  volume
245            z0ai  (:,:,jl)  = a_i  (:,:,  jl) * e1e2t(:,:)  ! Ice area
246            z0smi (:,:,jl)  = smv_i(:,:,  jl) * e1e2t(:,:)  ! Salt content
247            z0oi (:,:,jl)   = oa_i (:,:,  jl) * e1e2t(:,:)  ! Age content
248            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e1e2t(:,:)  ! Snow heat content
249            DO jk = 1, nlay_i
250               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
251            END DO
252            ! MV MP 2016
253            IF ( nn_pnd_scheme > 0 ) THEN
254               z0ap  (:,:,jl)  = a_ip (:,:,jl) * e1e2t(:,:) ! Melt pond fraction
255               z0vp  (:,:,jl)  = v_ip (:,:,jl) * e1e2t(:,:) ! Melt pond volume
256            ENDIF
257            ! END MV MP 2016
258
259         END DO
260
261
262         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
263            DO jt = 1, initad
264               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
265                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
266               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
267                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
268               DO jl = 1, jpl
269                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
270                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
271                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
272                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
273                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
274                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
275                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
276                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
277                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
278                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
279                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
280                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
281                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
282                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
283                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
284                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
285                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
286                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
287                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
288                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
289                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
290                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
291                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
292                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
293                  DO jk = 1, nlay_i                                                                !--- ice heat contents ---
294                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
295                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
296                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
297                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
298                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
299                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
300                  END DO
301                  ! MV MP 2016
302                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction --
303                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
304                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
305                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
306                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   --
307                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
308                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
309                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
310                  ! END MV MP 2016
311               END DO
312            END DO
313         ELSE
314            DO jt = 1, initad
315               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
316                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
317               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
318                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
319               DO jl = 1, jpl
320                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
321                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
322                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
323                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
324                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
325                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
326                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
327                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
328                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
329                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
330                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
331                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
332                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
333                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
334                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
335                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
336                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
337                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
338                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
339                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
340                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
341                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
342                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
343                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
344                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
345                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
346                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
347                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
348                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
349                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
350                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
351                  END DO
352                  ! MV MP 2016
353                  IF ( nn_pnd_scheme > 0 ) THEN
354                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction ---
355                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
356                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &
357                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
358                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   ---
359                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
360                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &
361                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
362                  ENDIF
363                  ! END MV MP 2016
364               END DO
365            END DO
366         ENDIF
367
368         !-------------------------------------------
369         ! Recover the properties from their contents
370         !-------------------------------------------
371         ato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:)
372         DO jl = 1, jpl
373            v_i  (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:)
374            v_s  (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:)
375            smv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:)
376            oa_i (:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:)
377            a_i  (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:)
378            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:)
379            DO jk = 1, nlay_i
380               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:)
381            END DO
382            ! MV MP 2016
383            IF ( nn_pnd_scheme > 0 ) THEN
384               a_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:)
385               v_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:)
386            ENDIF
387            ! END MV MP 2016
388         END DO
389
390         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
391         DO jl = 2, jpl
392            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
393         END DO
394         
395         CALL wrk_dealloc( jpi,jpj,            zarea )
396         CALL wrk_dealloc( jpi,jpj,1,          z0opw )
397         CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
398         ! MV MP 2016
399         CALL wrk_dealloc( jpi,jpj,jpl,        z0ap , z0vp )
400         ! END MV MP 2016
401         CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
402
403      END SELECT
404     
405      ! --- diags ---
406      DO jj = 1, jpj
407         DO ji = 1, jpi
408            diag_trp_ei (ji,jj) = ( SUM( e_i  (ji,jj,1:nlay_i,:) ) -  zeiold(ji,jj) ) * r1_rdtice
409            diag_trp_es (ji,jj) = ( SUM( e_s  (ji,jj,1:nlay_s,:) ) -  zesold(ji,jj) ) * r1_rdtice
410            diag_trp_smv(ji,jj) = ( SUM( smv_i(ji,jj,:)          ) - zsmvold(ji,jj) ) * r1_rdtice
411            diag_trp_vi (ji,jj) =   SUM(   v_i(ji,jj,:)            -  zviold(ji,jj,:) ) * r1_rdtice
412            diag_trp_vs (ji,jj) =   SUM(   v_s(ji,jj,:)            -  zvsold(ji,jj,:) ) * r1_rdtice
413         END DO
414      END DO
415     
416      IF( nn_limdyn == 2) THEN
417
418         ! zap small areas
419         CALL lim_var_zapsmall
420           
421         !--- Thickness correction in case too high --- !
422         DO jl = 1, jpl
423            DO jj = 1, jpj
424               DO ji = 1, jpi
425                 
426                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
427                     
428                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
429                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
430                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
431                     
432                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 
433                     
434                     IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
435                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
436                       
437                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
438                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
439                       
440                        ! small correction due to *rswitch for a_i
441                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl)
442                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl)
443                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl)
444                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
445                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
446
447                        ! MV MP 2016
448                        IF ( nn_pnd_scheme > 0 ) THEN
449                           a_ip (ji,jj,jl)        = rswitch * a_ip (ji,jj,jl)
450                           v_ip (ji,jj,jl)        = rswitch * v_ip (ji,jj,jl)
451                        ENDIF
452                        ! END MV MP 2016
453                                               
454                     ENDIF
455                     
456                  ENDIF
457               
458               END DO
459            END DO
460         END DO
461         ! -------------------------------------------------
462         
463         ! Force the upper limit of ht_i to always be < hi_max (99 m).
464         DO jj = 1, jpj
465            DO ji = 1, jpi
466               ! MV MP 2016
467               za_old = a_i(ji,jj,jpl)
468               ! END MV MP 2016
469               rswitch         = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
470               ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
471               a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
472               ! MV MP 2016
473               IF ( nn_pnd_scheme > 0 ) THEN
474                  ! correct pond fraction to avoid a_ip > a_i
475                  a_ip(ji,jj,jpl) = a_ip(ji,jj,jpl) * a_i(ji,jj,jpl) / MAX( za_old , epsi20 ) * rswitch
476               ENDIF
477               ! END MP 2016
478            END DO
479         END DO
480
481
482      ENDIF
483         
484      !------------------------------------------------------------
485      ! Impose a_i < amax if no ridging/rafting or in mono-category
486      !------------------------------------------------------------
487      !
488      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
489      IF ( nn_limdyn == 1 .OR. ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) ) THEN ! simple conservative piling, comparable with LIM2
490         DO jl = 1, jpl
491            DO jj = 1, jpj
492               DO ji = 1, jpi
493                  rswitch       = MAX( 0._wp, SIGN( 1._wp, at_i(ji,jj) - epsi20 ) )
494                  zda           = rswitch * MIN( rn_amax_2d(ji,jj) - at_i(ji,jj), 0._wp )  &
495                     &                    * a_i(ji,jj,jl) / MAX( at_i(ji,jj), epsi20 )
496                  a_i(ji,jj,jl) = a_i(ji,jj,jl) + zda
497               END DO
498            END DO
499         END DO
500      ENDIF
501     
502      ! --- agglomerate variables -----------------
503      vt_i(:,:) = SUM( v_i(:,:,:), dim=3 )
504      vt_s(:,:) = SUM( v_s(:,:,:), dim=3 )
505      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
506
507      ! MV MP 2016 (remove once we get rid of a_i_frac and ht_i)
508      IF ( nn_pnd_scheme > 0 ) THEN
509          at_ip(:,:) = SUM( a_ip(:,:,:), dim = 3 )
510          vt_ip(:,:) = SUM( v_ip(:,:,:), dim = 3 )
511      ENDIF
512      ! END MP 2016
513     
514      ! --- open water = 1 if at_i=0 --------------------------------
515      WHERE( at_i == 0._wp ) ato_i = 1._wp 
516     
517      ! conservation test
518      IF( ln_limdiachk )   CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
519       
520      ! -------------------------------------------------
521      ! control prints
522      ! -------------------------------------------------
523      IF( ln_limctl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
524      !
525      CALL wrk_dealloc( jpi,jpj,                            zatold, zeiold, zesold, zsmvold )
526      CALL wrk_dealloc( jpi,jpj,jpl,                        zhimax, zviold, zvsold )
527      !
528      IF( nn_timing == 1 )  CALL timing_stop('limtrp')
529      !
530   END SUBROUTINE lim_trp
531
532#else
533   !!----------------------------------------------------------------------
534   !!   Default option         Empty Module                No sea-ice model
535   !!----------------------------------------------------------------------
536CONTAINS
537   SUBROUTINE lim_trp        ! Empty routine
538   END SUBROUTINE lim_trp
539#endif
540   !!======================================================================
541END MODULE limtrp
542
Note: See TracBrowser for help on using the repository browser.