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

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

changes in style - part1 - (now the code looks better txs to Gurvan's comments)

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