source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90 @ 8426

Last change on this file since 8426 was 8426, checked in by clem, 3 years ago

last routine names to be changed

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