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_2.F90 in branches/2011/dev_r2802_LOCEAN10_agrif_lim/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2011/dev_r2802_LOCEAN10_agrif_lim/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90 @ 2804

Last change on this file since 2804 was 2804, checked in by rblod, 13 years ago

dev_r2802_LOCEAN10_agrif_lim: first implementation see ticket #848

  • Property svn:keywords set to Id
File size: 17.6 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
30# if defined key_agrif
31   USE agrif_lim2_interp ! nesting
32# endif
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   lim_trp_2   ! called by sbc_ice_lim_2
38
39   REAL(wp), PUBLIC ::   bound  = 0.e0          !: boundary condit. (0.0 no-slip, 1.0 free-slip)
40
41   REAL(wp)  ::   epsi06 = 1.e-06   ! constant values
42   REAL(wp)  ::   epsi03 = 1.e-03 
43   REAL(wp)  ::   epsi16 = 1.e-16 
44   REAL(wp)  ::   rzero  = 0.e0   
45   REAL(wp)  ::   rone   = 1.e0
46
47   !! * Substitution
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54
55CONTAINS
56
57   SUBROUTINE lim_trp_2( kt )
58      !!-------------------------------------------------------------------
59      !!                   ***  ROUTINE lim_trp_2 ***
60      !!                   
61      !! ** purpose : advection/diffusion process of sea ice
62      !!
63      !! ** method  : variables included in the process are scalar,   
64      !!     other values are considered as second order.
65      !!     For advection, a second order Prather scheme is used. 
66      !!
67      !! ** action :
68      !!---------------------------------------------------------------------
69      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
70      USE wrk_nemo, ONLY: zui_u  => wrk_2d_1, zvi_v => wrk_2d_2
71      USE wrk_nemo, ONLY: zsmice => wrk_2d_3, zsmsn => wrk_2d_4, zsma => wrk_2d_5
72      USE wrk_nemo, ONLY: zsmc0 => wrk_2d_6,  zsmc1 => wrk_2d_7, zsmc2 => wrk_2d_8, &
73                          zsmst => wrk_2d_9
74      !!
75      INTEGER, INTENT(in) ::   kt     ! number of iteration
76      !!
77      INTEGER  ::   ji, jj, jk   ! dummy loop indices
78      INTEGER  ::   initad       ! number of sub-timestep for the advection
79      REAL(wp) ::   zindb  , zindsn , zindic, zacrith   ! local scalars
80      REAL(wp) ::   zusvosn, zusvoic, zignm , zindhe    !   -      -
81      REAL(wp) ::   zvbord , zcfl   , zusnit            !   -      -
82      REAL(wp) ::   zrtt   , ztsn   , ztic1 , ztic2     !   -      -
83      !---------------------------------------------------------------------
84
85      IF( wrk_in_use(2, 1,2,3,4,5,6,7,8,9) ) THEN
86         CALL ctl_stop('lim_trp_2 : requested workspace arrays unavailable')   ;   RETURN
87      ENDIF
88
89      IF( kt == nit000  )   CALL lim_trp_init_2      ! Initialization (first time-step only)
90
91      zsmice(:,:) = area(:,:)
92      zsmsn (:,:) = area(:,:)
93      zsma  (:,:) = area(:,:)
94      zsmc0 (:,:) = area(:,:)
95      zsmc1 (:,:) = area(:,:)
96      zsmc2 (:,:) = area(:,:)
97      zsmst (:,:) = area(:,:)
98     
99      IF( ln_limdyn ) THEN
100         !-------------------------------------!
101         !   Advection of sea ice properties   !
102         !-------------------------------------!
103
104         ! ice velocities at ocean U- and V-points (zui_u,zvi_v)
105         ! ---------------------------------------
106         IF( lk_lim2_vp ) THEN      ! VP rheology : B-grid sea-ice dynamics (I-point ice velocity)
107            zvbord = 1._wp + ( 1._wp - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions       
108            DO jj = 1, jpjm1
109               DO ji = 1, jpim1   ! NO vector opt.
110                  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 ) )
111                  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 ) )
112               END DO
113            END DO
114            CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )      ! Lateral boundary conditions
115            !
116         ELSE                       ! EVP rheology : C-grid sea-ice dynamics (u- & v-points ice velocity)
117            zui_u(:,:) = u_ice(:,:)      ! EVP rheology: ice (u,v) at u- and v-points
118            zvi_v(:,:) = v_ice(:,:)
119         ENDIF
120
121         ! CFL test for stability
122         ! ----------------------
123         zcfl  = 0._wp
124         zcfl  = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, :     ) ) * rdt_ice / e1u(1:jpim1, :     ) ) )
125         zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) )
126         !
127         IF(lk_mpp)   CALL mpp_max( zcfl )
128         !
129         IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ', zcfl
130
131         ! content of properties
132         ! ---------------------
133         s0sn (:,:) =  hsnm(:,:)              * area  (:,:)  ! Snow volume.
134         s0ice(:,:) =  hicm(:,:)              * area  (:,:)  ! Ice volume.
135         s0a  (:,:) =  ( 1.0 - frld(:,:) )    * area  (:,:)  ! Surface covered by ice.
136         s0c0 (:,:) =  tbif(:,:,1) / rt0_snow * s0sn (:,:)  ! Heat content of the snow layer.
137         s0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * s0ice(:,:)  ! Heat content of the first ice layer.
138         s0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * s0ice(:,:)  ! Heat content of the second ice layer.
139         s0st (:,:) =  qstoif(:,:) / xlic     * s0a  (:,:)  ! Heat reservoir for brine pockets.
140         
141 
142         ! Advection (Prather scheme)
143         ! ---------
144         initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )   ! If ice drift field is too fast,         
145         zusnit = 1.0 / REAL( initad )                              ! split the ice time step in two
146         !
147         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN        !==  odd ice time step:  adv_x then adv_y  ==!
148            DO jk = 1, initad
149               CALL lim_adv_x_2( zusnit, zui_u, rone , zsmice, s0ice, sxice, sxxice, syice, syyice, sxyice )
150               CALL lim_adv_x_2( zusnit, zui_u, rone , zsmsn , s0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
151               CALL lim_adv_x_2( zusnit, zui_u, rone , zsma  , s0a  , sxa  , sxxa  , sya  , syya  , sxya   )
152               CALL lim_adv_x_2( zusnit, zui_u, rone , zsmc0 , s0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
153               CALL lim_adv_x_2( zusnit, zui_u, rone , zsmc1 , s0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
154               CALL lim_adv_x_2( zusnit, zui_u, rone , zsmc2 , s0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
155               CALL lim_adv_x_2( zusnit, zui_u, rone , zsmst , s0st , sxst , sxxst , syst , syyst , sxyst  )
156# if defined key_agrif
157               CALL Agrif_sadv_lim2(kt)
158# endif
159               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsmice, s0ice, sxice, sxxice, syice, syyice, sxyice )
160               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsmsn , s0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
161               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsma  , s0a  , sxa  , sxxa  , sya  , syya  , sxya   )
162               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsmc0 , s0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
163               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsmc1 , s0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
164               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsmc2 , s0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
165               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsmst , s0st , sxst , sxxst , syst , syyst , sxyst  )
166# if defined key_agrif
167               CALL Agrif_sadv_lim2(kt)
168# endif
169            END DO
170         ELSE                                                 !==  even ice time step:  adv_x then adv_y  ==!
171            DO jk = 1, initad
172               CALL lim_adv_y_2( zusnit, zvi_v, rone, zsmice, s0ice, sxice, sxxice, syice, syyice, sxyice )
173               CALL lim_adv_y_2( zusnit, zvi_v, rone, zsmsn , s0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
174               CALL lim_adv_y_2( zusnit, zvi_v, rone, zsma  , s0a  , sxa  , sxxa  , sya  , syya  , sxya   )
175               CALL lim_adv_y_2( zusnit, zvi_v, rone, zsmc0 , s0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
176               CALL lim_adv_y_2( zusnit, zvi_v, rone, zsmc1 , s0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
177               CALL lim_adv_y_2( zusnit, zvi_v, rone, zsmc2 , s0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
178               CALL lim_adv_y_2( zusnit, zvi_v, rone, zsmst , s0st , sxst , sxxst , syst , syyst , sxyst  )
179# if defined key_agrif
180               CALL Agrif_sadv_lim2(kt)
181# endif
182               CALL lim_adv_x_2( zusnit, zui_u, rzero , zsmice, s0ice, sxice, sxxice, syice, syyice, sxyice )
183               CALL lim_adv_x_2( zusnit, zui_u, rzero , zsmsn , s0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
184               CALL lim_adv_x_2( zusnit, zui_u, rzero , zsma  , s0a  , sxa  , sxxa  , sya  , syya  , sxya   )
185               CALL lim_adv_x_2( zusnit, zui_u, rzero , zsmc0 , s0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
186               CALL lim_adv_x_2( zusnit, zui_u, rzero , zsmc1 , s0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
187               CALL lim_adv_x_2( zusnit, zui_u, rzero , zsmc2 , s0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
188               CALL lim_adv_x_2( zusnit, zui_u, rzero , zsmst , s0st , sxst , sxxst , syst , syyst , sxyst  )
189# if defined key_agrif
190               CALL Agrif_sadv_lim2(kt)
191# endif
192            END DO
193         ENDIF
194                       
195         ! recover the properties from their contents
196         ! ------------------------------------------
197!!gm Define in limmsh one for all area = 1 /area  (CPU time saved !)
198         s0ice(:,:) = s0ice(:,:) / area(:,:)
199         s0sn (:,:) = s0sn (:,:) / area(:,:)
200         s0a  (:,:) = s0a  (:,:) / area(:,:)
201         s0c0 (:,:) = s0c0 (:,:) / area(:,:)
202         s0c1 (:,:) = s0c1 (:,:) / area(:,:)
203         s0c2 (:,:) = s0c2 (:,:) / area(:,:)
204         s0st (:,:) = s0st (:,:) / area(:,:)
205
206
207         !-------------------------------------!
208         !   Diffusion of sea ice properties   !
209         !-------------------------------------!
210
211         ! Masked eddy diffusivity coefficient at ocean U- and V-points
212         ! ------------------------------------------------------------
213         DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
214            DO ji = 1 , fs_jpim1   ! vector opt.
215               pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -s0a(ji  ,jj) ) ) )   &
216                  &        * ( 1.0 - MAX( rzero, SIGN( rone, -s0a(ji+1,jj) ) ) ) * ahiu(ji,jj)
217               pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -s0a(ji,jj  ) ) ) )   &
218                  &        * ( 1.0 - MAX( rzero, SIGN( rone,- s0a(ji,jj+1) ) ) ) * ahiv(ji,jj)
219            END DO
220         END DO
221!!gm more readable coding: (and avoid an error in F90 with sign of zero)
222!        DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
223!           DO ji = 1 , fs_jpim1   ! vector opt.
224!              IF( MIN( s0a(ji,jj) , s0a(ji+1,jj) ) == 0.e0 )   pahu(ji,jj) = 0._wp
225!              IF( MIN( s0a(ji,jj) , s0a(ji,jj+1) ) == 0.e0 )   pahv(ji,jj) = 0._wp
226!           END DO
227!        END DO
228!!gm end
229
230         ! diffusion
231         ! ---------
232         CALL lim_hdf_2( s0ice )
233         CALL lim_hdf_2( s0sn  )
234         CALL lim_hdf_2( s0a   )
235         CALL lim_hdf_2( s0c0  )
236         CALL lim_hdf_2( s0c1  )
237         CALL lim_hdf_2( s0c2  )
238         CALL lim_hdf_2( s0st  )
239
240!!gm see comment this can be skipped
241         s0ice(:,:) = MAX( rzero, s0ice(:,:) * area(:,:) )    !!bug:  useless
242         s0sn (:,:) = MAX( rzero, s0sn (:,:) * area(:,:) )    !!bug:  cf /area  just below
243         s0a  (:,:) = MAX( rzero, s0a  (:,:) * area(:,:) )    !! caution: the suppression of the 2 changes
244         s0c0 (:,:) = MAX( rzero, s0c0 (:,:) * area(:,:) )    !! the last digit of the results
245         s0c1 (:,:) = MAX( rzero, s0c1 (:,:) * area(:,:) )
246         s0c2 (:,:) = MAX( rzero, s0c2 (:,:) * area(:,:) )
247         s0st (:,:) = MAX( rzero, s0st (:,:) * area(:,:) )
248
249
250         !-------------------------------------------------------------------!
251         !   Updating and limitation of sea ice properties after transport   !
252         !-------------------------------------------------------------------!
253         DO jj = 1, jpj
254            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH
255            DO ji = 1, jpi
256               !
257               ! Recover mean values over the grid squares.
258               s0sn (ji,jj) = MAX( rzero, s0sn (ji,jj)/area(ji,jj) )
259               s0ice(ji,jj) = MAX( rzero, s0ice(ji,jj)/area(ji,jj) )
260               s0a  (ji,jj) = MAX( rzero, s0a  (ji,jj)/area(ji,jj) )
261               s0c0 (ji,jj) = MAX( rzero, s0c0 (ji,jj)/area(ji,jj) )
262               s0c1 (ji,jj) = MAX( rzero, s0c1 (ji,jj)/area(ji,jj) )
263               s0c2 (ji,jj) = MAX( rzero, s0c2 (ji,jj)/area(ji,jj) )
264               s0st (ji,jj) = MAX( rzero, s0st (ji,jj)/area(ji,jj) )
265
266               ! Recover in situ values.
267               zindb         = MAX( rzero, SIGN( rone, s0a(ji,jj) - epsi06 ) )
268               zacrith       = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) )
269               s0a (ji,jj)  = zindb * MIN( s0a(ji,jj), zacrith )
270               hsnif(ji,jj)  = zindb * ( s0sn(ji,jj) /MAX( s0a(ji,jj), epsi16 ) )
271               hicif(ji,jj)  = zindb * ( s0ice(ji,jj)/MAX( s0a(ji,jj), epsi16 ) )
272               zindsn        = MAX( rzero, SIGN( rone, hsnif(ji,jj) - epsi06 ) )
273               zindic        = MAX( rzero, SIGN( rone, hicif(ji,jj) - epsi03 ) )
274               zindb         = MAX( zindsn, zindic )
275               s0a (ji,jj)  = zindb * s0a(ji,jj)
276               frld (ji,jj)  = 1.0 - s0a(ji,jj)
277               hsnif(ji,jj)  = zindsn * hsnif(ji,jj)
278               hicif(ji,jj)  = zindic * hicif(ji,jj)
279               zusvosn       = 1.0/MAX( hsnif(ji,jj) * s0a(ji,jj), epsi16 )
280               zusvoic       = 1.0/MAX( hicif(ji,jj) * s0a(ji,jj), epsi16 )
281               zignm         = MAX( rzero,  SIGN( rone, hsndif - hsnif(ji,jj) ) )
282               zrtt          = 173.15 * rone 
283               ztsn          =          zignm   * tbif(ji,jj,1)  &
284                              + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * s0c0(ji,jj)) , tfu(ji,jj) ) 
285               ztic1          = MIN( MAX( zrtt, rt0_ice * zusvoic * s0c1(ji,jj) ) , tfu(ji,jj) )
286               ztic2          = MIN( MAX( zrtt, rt0_ice * zusvoic * s0c2(ji,jj) ) , tfu(ji,jj) )
287 
288               tbif(ji,jj,1) = zindsn * ztsn  + ( 1.0 - zindsn ) * tfu(ji,jj)               
289               tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj)
290               tbif(ji,jj,3) = zindic * ztic2 + ( 1.0 - zindic ) * tfu(ji,jj)
291               qstoif(ji,jj) = zindb  * xlic * s0st(ji,jj) /  MAX( s0a(ji,jj), epsi16 )
292            END DO
293         END DO
294         !
295      ENDIF
296      !
297# if defined key_agrif
298!RB : less precise but could be enough
299!     and avoid above calls to Agrif_sadv_lim2
300              CALL Agrif_adv_lim2(kt)
301# endif
302      !
303      IF( wrk_not_released(2, 1,2,3,4,5,6,7,8,9) )   &
304         &   CALL ctl_stop('lim_trp_2 : failed to release workspace arrays')
305      !
306   END SUBROUTINE lim_trp_2
307
308
309   SUBROUTINE lim_trp_init_2
310      !!-------------------------------------------------------------------
311      !!                  ***  ROUTINE lim_trp_init_2  ***
312      !!
313      !! ** Purpose :   initialization of ice advection parameters
314      !!
315      !! ** Method  :   Read the namicetrp namelist and check the parameter
316      !!              values called at the first timestep (nit000)
317      !!
318      !! ** input   :   Namelist namicetrp
319      !!-------------------------------------------------------------------
320      NAMELIST/namicetrp/ bound
321      !!-------------------------------------------------------------------
322      !
323      REWIND ( numnam_ice )      ! Read Namelist namicetrp
324      READ   ( numnam_ice  , namicetrp )
325      IF(lwp) THEN
326         WRITE(numout,*)
327         WRITE(numout,*) 'lim_trp_init_2 : Ice parameters for advection '
328         WRITE(numout,*) '~~~~~~~~~~~~~~'
329         WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound
330      ENDIF
331      !
332   END SUBROUTINE lim_trp_init_2
333
334#else
335   !!----------------------------------------------------------------------
336   !!   Default option         Empty Module                No sea-ice model
337   !!----------------------------------------------------------------------
338CONTAINS
339   SUBROUTINE lim_trp_2        ! Empty routine
340   END SUBROUTINE lim_trp_2
341#endif
342
343   !!======================================================================
344END MODULE limtrp_2
Note: See TracBrowser for help on using the repository browser.