source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

File size: 17.7 KB
Line 
1MODULE limtrp_2
2   !!======================================================================
3   !!                       ***  MODULE limtrp_2   ***
4   !! LIM 2.0 transport ice model : sea-ice advection/diffusion
5   !!======================================================================
6   !! History :  LIM  !  2000-01 (UCL)  Original code
7   !!            2.0  !  2001-05 (G. Madec, R. Hordoir) opa norm
8   !!             -   !  2004-01 (G. Madec, C. Ethe)  F90, mpp
9   !!            3.3  !  2009-05  (G. Garric, C. Bricaud) addition of the lim2_evp case
10   !!----------------------------------------------------------------------
11#if defined key_lim2
12   !!----------------------------------------------------------------------
13   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_trp_2      : advection/diffusion process of sea ice
16   !!   lim_trp_init_2 : initialization and namelist read
17   !!----------------------------------------------------------------------
18   USE phycst          ! physical constant
19   USE sbc_oce         ! ocean surface boundary condition
20   USE dom_oce         ! ocean domain
21   USE in_out_manager  ! I/O manager
22   USE dom_ice_2       ! LIM-2 domain
23   USE ice_2           ! LIM-2 variables
24   USE limistate_2     ! LIM-2 initial state
25   USE limadv_2        ! LIM-2 advection
26   USE limhdf_2        ! LIM-2 horizontal diffusion
27   USE lbclnk          ! lateral boundary conditions -- MPP exchanges
28   USE lib_mpp         ! MPP library
29   USE wrk_nemo        ! work arrays
30# if defined key_agrif
31   USE agrif_lim2_interp ! nesting
32# endif
33   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   lim_trp_2   ! called by sbc_ice_lim_2
39
40   REAL(wp), PUBLIC ::   bound      !: boundary condit. (0.0 no-slip, 1.0 free-slip)
41
42   REAL(wp)  ::   epsi06 = 1.e-06   ! constant values
43   REAL(wp)  ::   epsi03 = 1.e-03 
44   REAL(wp)  ::   epsi16 = 1.e-16 
45   REAL(wp)  ::   rzero  = 0.e0   
46   REAL(wp)  ::   rone   = 1.e0
47
48   !! * Substitution
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
52   !! $Id$
53   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55
56CONTAINS
57
58   SUBROUTINE lim_trp_2( kt )
59      !!-------------------------------------------------------------------
60      !!                   ***  ROUTINE lim_trp_2 ***
61      !!                   
62      !! ** purpose : advection/diffusion process of sea ice
63      !!
64      !! ** method  : variables included in the process are scalar,   
65      !!     other values are considered as second order.
66      !!     For advection, a second order Prather scheme is used. 
67      !!
68      !! ** action :
69      !!---------------------------------------------------------------------
70      INTEGER, INTENT(in) ::   kt     ! number of iteration
71      !!
72      INTEGER  ::   ji, jj, jk   ! dummy loop indices
73      INTEGER  ::   initad       ! number of sub-timestep for the advection
74      REAL(wp) ::   zindb  , zindsn , zindic, zacrith   ! local scalars
75      REAL(wp) ::   zusvosn, zusvoic, zignm , zindhe    !   -      -
76      REAL(wp) ::   zvbord , zcfl   , zusnit            !   -      -
77      REAL(wp) ::   zrtt   , ztsn   , ztic1 , ztic2     !   -      -
78      REAL(wp), POINTER, DIMENSION(:,:)  ::   zui_u , zvi_v , zsm             ! 2D workspace
79      REAL(wp), POINTER, DIMENSION(:,:)  ::   zs0ice, zs0sn , zs0a            !  -      -
80      REAL(wp), POINTER, DIMENSION(:,:)  ::   zs0c0 , zs0c1 , zs0c2 , zs0st   !  -      -
81      !---------------------------------------------------------------------
82
83      CALL wrk_alloc( jpi, jpj, zui_u , zvi_v , zsm, zs0ice, zs0sn , zs0a, zs0c0 , zs0c1 , zs0c2 , zs0st )
84
85      IF( kt == nit000  )   CALL lim_trp_init_2      ! Initialization (first time-step only)
86
87# if defined key_agrif
88      CALL agrif_trp_lim2_load      ! First interpolation
89# endif
90
91      zsm(:,:) = area(:,:)
92     
93      IF( ln_limdyn ) THEN
94         !-------------------------------------!
95         !   Advection of sea ice properties   !
96         !-------------------------------------!
97
98         ! ice velocities at ocean U- and V-points (zui_u,zvi_v)
99         ! ---------------------------------------
100         IF( lk_lim2_vp ) THEN      ! VP rheology : B-grid sea-ice dynamics (I-point ice velocity)
101            zvbord = 1._wp + ( 1._wp - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions       
102            DO jj = 1, jpjm1
103               DO ji = 1, jpim1   ! NO vector opt.
104                  zui_u(ji,jj) = ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj)+tmu(ji+1,jj+1), zvbord ) )
105                  zvi_v(ji,jj) = ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji,jj+1)+tmu(ji+1,jj+1), zvbord ) )
106               END DO
107            END DO
108            CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )      ! Lateral boundary conditions
109            !
110         ELSE                       ! EVP rheology : C-grid sea-ice dynamics (u- & v-points ice velocity)
111            zui_u(:,:) = u_ice(:,:)      ! EVP rheology: ice (u,v) at u- and v-points
112            zvi_v(:,:) = v_ice(:,:)
113         ENDIF
114
115         ! CFL test for stability
116         ! ----------------------
117         zcfl  = 0._wp
118         zcfl  = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, :     ) ) * rdt_ice / e1u(1:jpim1, :     ) ) )
119         zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) )
120         !
121         IF(lk_mpp)   CALL mpp_max( zcfl )
122         !
123         IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ', zcfl
124
125         ! content of properties
126         ! ---------------------
127         zs0sn (:,:) =  hsnm(:,:)              * area  (:,:)  ! Snow volume.
128         zs0ice(:,:) =  hicm(:,:)              * area  (:,:)  ! Ice volume.
129         zs0a  (:,:) =  ( 1.0 - frld(:,:) )    * area  (:,:)  ! Surface covered by ice.
130         zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn (:,:)  ! Heat content of the snow layer.
131         zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)  ! Heat content of the first ice layer.
132         zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)  ! Heat content of the second ice layer.
133         zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a  (:,:)  ! Heat reservoir for brine pockets.
134         
135 
136         ! Advection (Prather scheme)
137         ! ---------
138         initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )   ! If ice drift field is too fast,         
139         zusnit = 1.0 / REAL( initad )                              ! split the ice time step in two
140         !
141         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN        !==  odd ice time step:  adv_x then adv_y  ==!
142            DO jk = 1, initad
143               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
144               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
145               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
146               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
147               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
148               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
149               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
150               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
151               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
152               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
153               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
154               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
155               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
156               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
157            END DO
158         ELSE                                                 !==  even ice time step:  adv_x then adv_y  ==!
159            DO jk = 1, initad
160               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
161               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
162               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
163               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
164               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
165               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
166               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
167               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
168               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
169               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
170               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
171               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
172               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
173               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
174            END DO
175         ENDIF
176                       
177         ! recover the properties from their contents
178         ! ------------------------------------------
179!!gm Define in limmsh one for all area = 1 /area  (CPU time saved !)
180         zs0ice(:,:) = zs0ice(:,:) / area(:,:)
181         zs0sn (:,:) = zs0sn (:,:) / area(:,:)
182         zs0a  (:,:) = zs0a  (:,:) / area(:,:)
183         zs0c0 (:,:) = zs0c0 (:,:) / area(:,:)
184         zs0c1 (:,:) = zs0c1 (:,:) / area(:,:)
185         zs0c2 (:,:) = zs0c2 (:,:) / area(:,:)
186         zs0st (:,:) = zs0st (:,:) / area(:,:)
187
188
189         !-------------------------------------!
190         !   Diffusion of sea ice properties   !
191         !-------------------------------------!
192
193         ! Masked eddy diffusivity coefficient at ocean U- and V-points
194         ! ------------------------------------------------------------
195         DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
196            DO ji = 1 , fs_jpim1   ! vector opt.
197               pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji  ,jj) ) ) )   &
198                  &        * ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj) ) ) ) * ahiu(ji,jj)
199               pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji,jj  ) ) ) )   &
200                  &        * ( 1.0 - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1) ) ) ) * ahiv(ji,jj)
201            END DO
202         END DO
203!!gm more readable coding: (and avoid an error in F90 with sign of zero)
204!        DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
205!           DO ji = 1 , fs_jpim1   ! vector opt.
206!              IF( MIN( zs0a(ji,jj) , zs0a(ji+1,jj) ) == 0.e0 )   pahu(ji,jj) = 0._wp
207!              IF( MIN( zs0a(ji,jj) , zs0a(ji,jj+1) ) == 0.e0 )   pahv(ji,jj) = 0._wp
208!           END DO
209!        END DO
210!!gm end
211
212         ! diffusion
213         ! ---------
214         CALL lim_hdf_2( zs0ice )
215         CALL lim_hdf_2( zs0sn  )
216         CALL lim_hdf_2( zs0a   )
217         CALL lim_hdf_2( zs0c0  )
218         CALL lim_hdf_2( zs0c1  )
219         CALL lim_hdf_2( zs0c2  )
220         CALL lim_hdf_2( zs0st  )
221
222!!gm see comment this can be skipped
223         zs0ice(:,:) = MAX( rzero, zs0ice(:,:) * area(:,:) )    !!bug:  useless
224         zs0sn (:,:) = MAX( rzero, zs0sn (:,:) * area(:,:) )    !!bug:  cf /area  just below
225         zs0a  (:,:) = MAX( rzero, zs0a  (:,:) * area(:,:) )    !! caution: the suppression of the 2 changes
226         zs0c0 (:,:) = MAX( rzero, zs0c0 (:,:) * area(:,:) )    !! the last digit of the results
227         zs0c1 (:,:) = MAX( rzero, zs0c1 (:,:) * area(:,:) )
228         zs0c2 (:,:) = MAX( rzero, zs0c2 (:,:) * area(:,:) )
229         zs0st (:,:) = MAX( rzero, zs0st (:,:) * area(:,:) )
230
231
232         !-------------------------------------------------------------------!
233         !   Updating and limitation of sea ice properties after transport   !
234         !-------------------------------------------------------------------!
235         DO jj = 1, jpj
236            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH
237            DO ji = 1, jpi
238               !
239               ! Recover mean values over the grid squares.
240               zs0sn (ji,jj) = MAX( rzero, zs0sn (ji,jj)/area(ji,jj) )
241               zs0ice(ji,jj) = MAX( rzero, zs0ice(ji,jj)/area(ji,jj) )
242               zs0a  (ji,jj) = MAX( rzero, zs0a  (ji,jj)/area(ji,jj) )
243               zs0c0 (ji,jj) = MAX( rzero, zs0c0 (ji,jj)/area(ji,jj) )
244               zs0c1 (ji,jj) = MAX( rzero, zs0c1 (ji,jj)/area(ji,jj) )
245               zs0c2 (ji,jj) = MAX( rzero, zs0c2 (ji,jj)/area(ji,jj) )
246               zs0st (ji,jj) = MAX( rzero, zs0st (ji,jj)/area(ji,jj) )
247
248               ! Recover in situ values.
249               zindb         = MAX( rzero, SIGN( rone, zs0a(ji,jj) - epsi06 ) )
250               zacrith       = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) )
251               zs0a (ji,jj)  = zindb * MIN( zs0a(ji,jj), zacrith )
252               hsnif(ji,jj)  = zindb * ( zs0sn(ji,jj) /MAX( zs0a(ji,jj), epsi16 ) )
253               hicif(ji,jj)  = zindb * ( zs0ice(ji,jj)/MAX( zs0a(ji,jj), epsi16 ) )
254               zindsn        = MAX( rzero, SIGN( rone, hsnif(ji,jj) - epsi06 ) )
255               zindic        = MAX( rzero, SIGN( rone, hicif(ji,jj) - epsi03 ) )
256               zindb         = MAX( zindsn, zindic )
257               zs0a (ji,jj)  = zindb * zs0a(ji,jj)
258               frld (ji,jj)  = 1.0 - zs0a(ji,jj)
259               hsnif(ji,jj)  = zindsn * hsnif(ji,jj)
260               hicif(ji,jj)  = zindic * hicif(ji,jj)
261               zusvosn       = 1.0/MAX( hsnif(ji,jj) * zs0a(ji,jj), epsi16 )
262               zusvoic       = 1.0/MAX( hicif(ji,jj) * zs0a(ji,jj), epsi16 )
263               zignm         = MAX( rzero,  SIGN( rone, hsndif - hsnif(ji,jj) ) )
264               zrtt          = 173.15 * rone 
265               ztsn          =          zignm   * tbif(ji,jj,1)  &
266                              + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) ) 
267               ztic1          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c1(ji,jj) ) , tfu(ji,jj) )
268               ztic2          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c2(ji,jj) ) , tfu(ji,jj) )
269 
270               tbif(ji,jj,1) = zindsn * ztsn  + ( 1.0 - zindsn ) * tfu(ji,jj)               
271               tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj)
272               tbif(ji,jj,3) = zindic * ztic2 + ( 1.0 - zindic ) * tfu(ji,jj)
273               qstoif(ji,jj) = zindb  * xlic * zs0st(ji,jj) /  MAX( zs0a(ji,jj), epsi16 )
274            END DO
275         END DO
276         !
277      ENDIF
278      !
279# if defined key_agrif
280      CALL agrif_trp_lim2      ! Fill boundaries of the fine grid
281# endif
282      !
283      CALL wrk_dealloc( jpi, jpj, zui_u , zvi_v , zsm, zs0ice, zs0sn , zs0a, zs0c0 , zs0c1 , zs0c2 , zs0st )
284      !
285   END SUBROUTINE lim_trp_2
286
287
288   SUBROUTINE lim_trp_init_2
289      !!-------------------------------------------------------------------
290      !!                  ***  ROUTINE lim_trp_init_2  ***
291      !!
292      !! ** Purpose :   initialization of ice advection parameters
293      !!
294      !! ** Method  :   Read the namicetrp namelist and check the parameter
295      !!              values called at the first timestep (nit000)
296      !!
297      !! ** input   :   Namelist namicetrp
298      !!-------------------------------------------------------------------
299      INTEGER  ::   ios                 ! Local integer output status for namelist read
300      NAMELIST/namicetrp/ bound
301      !!-------------------------------------------------------------------
302                   
303      REWIND( numnam_ice_ref )              ! Namelist namicetrp in reference namelist : Ice advection
304      READ  ( numnam_ice_ref, namicetrp, IOSTAT = ios, ERR = 901)
305901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicetrp in reference namelist', lwp )
306
307      REWIND( numnam_ice_cfg )              ! Namelist namicetrp in configuration namelist : Ice advection
308      READ  ( numnam_ice_cfg, namicetrp, IOSTAT = ios, ERR = 902 )
309902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicetrp in configuration namelist', lwp )
310      IF(lwm) WRITE ( numoni, namicetrp )
311
312      IF(lwp) THEN
313         WRITE(numout,*)
314         WRITE(numout,*) 'lim_trp_init_2 : Ice parameters for advection '
315         WRITE(numout,*) '~~~~~~~~~~~~~~'
316         WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound
317      ENDIF
318      !
319   END SUBROUTINE lim_trp_init_2
320
321#else
322   !!----------------------------------------------------------------------
323   !!   Default option         Empty Module                No sea-ice model
324   !!----------------------------------------------------------------------
325CONTAINS
326   SUBROUTINE lim_trp_2        ! Empty routine
327   END SUBROUTINE lim_trp_2
328#endif
329
330   !!======================================================================
331END MODULE limtrp_2
Note: See TracBrowser for help on using the repository browser.