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/devmercator2010/NEMO/LIM_SRC_2 – NEMO

source: branches/devmercator2010/NEMO/LIM_SRC_2/limtrp_2.F90 @ 2076

Last change on this file since 2076 was 2076, checked in by cbricaud, 14 years ago

add change dev_1784_EVP
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 15.5 KB
Line 
1MODULE limtrp_2
2   !!======================================================================
3   !!                       ***  MODULE limtrp_2   ***
4   !! LIM 2.0 transport ice model : sea-ice advection/diffusion
5   !!======================================================================
6#if defined key_lim2
7   !!----------------------------------------------------------------------
8   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
9   !!----------------------------------------------------------------------
10   !!   lim_trp_2      : advection/diffusion process of sea ice
11   !!   lim_trp_init_2 : initialization and namelist read
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE phycst
15   USE dom_oce
16   USE in_out_manager  ! I/O manager
17   USE dom_ice_2
18   USE ice_2
19   USE limistate_2
20   USE limadv_2
21   USE limhdf_2
22   USE lbclnk
23   USE lib_mpp
24
25   IMPLICIT NONE
26   PRIVATE
27
28   !! * Routine accessibility
29   PUBLIC lim_trp_2     ! called by sbc_ice_lim_2
30
31   !! * Shared module variables
32   REAL(wp), PUBLIC  ::   &  !:
33      bound  = 0.e0          !: boundary condit. (0.0 no-slip, 1.0 free-slip)
34
35   !! * Module variables
36   REAL(wp)  ::           &  ! constant values
37      epsi06 = 1.e-06  ,  &
38      epsi03 = 1.e-03  ,  &
39      epsi16 = 1.e-16  ,  &
40      rzero  = 0.e0    ,  &
41      rone   = 1.e0
42
43   !! * Substitution
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)
47   !! $Id$
48   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
49   !!----------------------------------------------------------------------
50
51CONTAINS
52
53   SUBROUTINE lim_trp_2( kt )
54      !!-------------------------------------------------------------------
55      !!                   ***  ROUTINE lim_trp_2 ***
56      !!                   
57      !! ** purpose : advection/diffusion process of sea ice
58      !!
59      !! ** method  : variables included in the process are scalar,   
60      !!     other values are considered as second order.
61      !!     For advection, a second order Prather scheme is used. 
62      !!
63      !! ** action :
64      !!
65      !! History :
66      !!   1.0  !  00-01 (LIM)  Original code
67      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
68      !!   2.0  !  04-01 (G. Madec, C. Ethe)  F90, mpp
69      !!   3.3  !  09-05  (G.Garric) addition of the lim2_evp case
70      !!---------------------------------------------------------------------
71      INTEGER, INTENT(in) ::   kt     ! number of iteration
72
73      INTEGER  ::   ji, jj, jk,   &  ! dummy loop indices
74         &          initad           ! number of sub-timestep for the advection
75
76      REAL(wp) ::  &                             
77         zindb  ,  &
78         zacrith, &
79         zindsn , &
80         zindic , &
81         zusvosn, &
82         zusvoic, &
83         zignm  , &
84         zindhe , &
85         zvbord , &
86         zcfl   , &
87         zusnit , &
88         zrtt, ztsn, ztic1, ztic2
89
90      REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace
91         zui_u , zvi_v , zsm   ,         &
92         zs0ice, zs0sn , zs0a  ,         &
93         zs0c0 , zs0c1 , zs0c2 ,         &
94         zs0st
95      !---------------------------------------------------------------------
96
97      IF( kt == nit000  )   CALL lim_trp_init_2      ! Initialization (first time-step only)
98
99      zsm(:,:) = area(:,:)
100     
101      IF( ln_limdyn ) THEN
102         !-------------------------------------!
103         !   Advection of sea ice properties   !
104         !-------------------------------------!
105
106         ! ice velocities at ocean U- and V-points (zui_u,zvi_v)
107         ! ---------------------------------------
108         ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions.       
109         zvbord = 1.0 + ( 1.0 - bound )
110#if defined key_lim2_vp
111         DO jj = 1, jpjm1
112            DO ji = 1, jpim1   ! NO vector opt.
113               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 ) )
114               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 ) )
115            END DO
116         END DO
117         ! Lateral boundary conditions on zui_u, zvi_v
118         CALL lbc_lnk( zui_u, 'U', -1. )
119         CALL lbc_lnk( zvi_v, 'V', -1. )
120#else
121        zui_u(:,:)=u_ice(:,:)
122        zvi_v(:,:)=v_ice(:,:)
123#endif
124
125         ! CFL test for stability
126         ! ----------------------
127         zcfl  = 0.e0
128         zcfl  = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, :     ) ) * rdt_ice / e1u(1:jpim1, :     ) ) )
129         zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) )
130
131         IF (lk_mpp ) CALL mpp_max(zcfl)
132
133         IF ( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl
134
135         ! content of properties
136         ! ---------------------
137         zs0sn (:,:) =  hsnm(:,:) * area(:,:)                 ! Snow volume.
138         zs0ice(:,:) =  hicm (:,:) * area(:,:)                ! Ice volume.
139         zs0a  (:,:) =  ( 1.0 - frld(:,:) ) * area(:,:)       ! Surface covered by ice.
140         zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn(:,:)   ! Heat content of the snow layer.
141         zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)  ! Heat content of the first ice layer.
142         zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)  ! Heat content of the second ice layer.
143         zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a(:,:)    ! Heat reservoir for brine pockets.
144         
145 
146         ! Advection
147         ! ---------
148         ! If ice drift field is too fast, use an appropriate time step for advection.         
149         initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )
150         zusnit = 1.0 / REAL( initad ) 
151         
152         IF ( MOD( nday , 2 ) == 0) THEN
153            DO jk = 1,initad
154               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
155               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
156               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
157               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
158               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
159               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
160               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
161               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
162               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
163               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
164               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
165               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
166               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
167               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
168            END DO
169         ELSE
170            DO jk = 1, initad
171               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
172               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
173               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
174               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
175               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
176               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
177               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
178               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
179               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
180               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
181               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
182               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
183               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
184               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
185            END DO
186         ENDIF
187                       
188         ! recover the properties from their contents
189         ! ------------------------------------------
190         zs0ice(:,:) = zs0ice(:,:) / area(:,:)
191         zs0sn (:,:) = zs0sn (:,:) / area(:,:)
192         zs0a  (:,:) = zs0a  (:,:) / area(:,:)
193         zs0c0 (:,:) = zs0c0 (:,:) / area(:,:)
194         zs0c1 (:,:) = zs0c1 (:,:) / area(:,:)
195         zs0c2 (:,:) = zs0c2 (:,:) / area(:,:)
196         zs0st (:,:) = zs0st (:,:) / area(:,:)
197
198
199         !-------------------------------------!
200         !   Diffusion of sea ice properties   !
201         !-------------------------------------!
202
203         ! Masked eddy diffusivity coefficient at ocean U- and V-points
204         ! ------------------------------------------------------------
205         DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
206            DO ji = 1 , fs_jpim1   ! vector opt.
207               pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji  ,jj) ) ) )   &
208                  &        * ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj) ) ) ) * ahiu(ji,jj)
209               pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji,jj  ) ) ) )   &
210                  &        * ( 1.0 - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1) ) ) ) * ahiv(ji,jj)
211            END DO
212         END DO
213
214         ! diffusion
215         ! ---------
216         CALL lim_hdf_2( zs0ice )
217         CALL lim_hdf_2( zs0sn  )
218         CALL lim_hdf_2( zs0a   )
219         CALL lim_hdf_2( zs0c0  )
220         CALL lim_hdf_2( zs0c1  )
221         CALL lim_hdf_2( zs0c2  )
222         CALL lim_hdf_2( zs0st  )
223
224         zs0ice(:,:) = MAX( rzero, zs0ice(:,:) * area(:,:) )    !!bug:  est-ce utile
225         zs0sn (:,:) = MAX( rzero, zs0sn (:,:) * area(:,:) )    !!bug:  cf /area  juste apres
226         zs0a  (:,:) = MAX( rzero, zs0a  (:,:) * area(:,:) )    !! suppression des 2 change le resultat...
227         zs0c0 (:,:) = MAX( rzero, zs0c0 (:,:) * area(:,:) )
228         zs0c1 (:,:) = MAX( rzero, zs0c1 (:,:) * area(:,:) )
229         zs0c2 (:,:) = MAX( rzero, zs0c2 (:,:) * area(:,:) )
230         zs0st (:,:) = MAX( rzero, zs0st (:,:) * area(:,:) )
231
232
233         ! -------------------------------------------------------------------!
234         !   Up-dating and limitation of sea ice properties after transport   !
235         ! -------------------------------------------------------------------!
236
237         ! Up-dating and limitation of sea ice properties after transport.
238         DO jj = 1, jpj
239!!!iii      zindhe = REAL( MAX( 0, isign(1, jj - njeq ) ) )              !ibug mpp  !!bugmpp  njeq!
240            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH
241            DO ji = 1, jpi
242
243               ! Recover mean values over the grid squares.
244               zs0sn (ji,jj) = MAX( rzero, zs0sn (ji,jj)/area(ji,jj) )
245               zs0ice(ji,jj) = MAX( rzero, zs0ice(ji,jj)/area(ji,jj) )
246               zs0a  (ji,jj) = MAX( rzero, zs0a  (ji,jj)/area(ji,jj) )
247               zs0c0 (ji,jj) = MAX( rzero, zs0c0 (ji,jj)/area(ji,jj) )
248               zs0c1 (ji,jj) = MAX( rzero, zs0c1 (ji,jj)/area(ji,jj) )
249               zs0c2 (ji,jj) = MAX( rzero, zs0c2 (ji,jj)/area(ji,jj) )
250               zs0st (ji,jj) = MAX( rzero, zs0st (ji,jj)/area(ji,jj) )
251
252               ! Recover in situ values.
253               zindb         = MAX( rzero, SIGN( rone, zs0a(ji,jj) - epsi06 ) )
254               zacrith       = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) )
255               zs0a (ji,jj)  = zindb * MIN( zs0a(ji,jj), zacrith )
256               hsnif(ji,jj)  = zindb * ( zs0sn(ji,jj) /MAX( zs0a(ji,jj), epsi16 ) )
257               hicif(ji,jj)  = zindb * ( zs0ice(ji,jj)/MAX( zs0a(ji,jj), epsi16 ) )
258               zindsn        = MAX( rzero, SIGN( rone, hsnif(ji,jj) - epsi06 ) )
259               zindic        = MAX( rzero, SIGN( rone, hicif(ji,jj) - epsi03 ) )
260               zindb         = MAX( zindsn, zindic )
261               zs0a (ji,jj)  = zindb * zs0a(ji,jj)
262               frld (ji,jj)  = 1.0 - zs0a(ji,jj)
263               hsnif(ji,jj)  = zindsn * hsnif(ji,jj)
264               hicif(ji,jj)  = zindic * hicif(ji,jj)
265               zusvosn       = 1.0/MAX( hsnif(ji,jj) * zs0a(ji,jj), epsi16 )
266               zusvoic       = 1.0/MAX( hicif(ji,jj) * zs0a(ji,jj), epsi16 )
267               zignm         = MAX( rzero,  SIGN( rone, hsndif - hsnif(ji,jj) ) )
268               zrtt          = 173.15 * rone 
269               ztsn          =          zignm   * tbif(ji,jj,1)  &
270                              + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) ) 
271               ztic1          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c1(ji,jj) ) , tfu(ji,jj) )
272               ztic2          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c2(ji,jj) ) , tfu(ji,jj) )
273 
274               tbif(ji,jj,1) = zindsn * ztsn  + ( 1.0 - zindsn ) * tfu(ji,jj)               
275               tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj)
276               tbif(ji,jj,3) = zindic * ztic2 + ( 1.0 - zindic ) * tfu(ji,jj)
277               qstoif(ji,jj) = zindb  * xlic * zs0st(ji,jj) /  MAX( zs0a(ji,jj), epsi16 )
278            END DO
279         END DO
280         
281      ENDIF
282     
283   END SUBROUTINE lim_trp_2
284
285
286   SUBROUTINE lim_trp_init_2
287      !!-------------------------------------------------------------------
288      !!                  ***  ROUTINE lim_trp_init_2  ***
289      !!
290      !! ** Purpose :   initialization of ice advection parameters
291      !!
292      !! ** Method  : Read the namicetrp namelist and check the parameter
293      !!       values called at the first timestep (nit000)
294      !!
295      !! ** input   :   Namelist namicetrp
296      !!
297      !! history :
298      !!   2.0  !  03-08 (C. Ethe)  Original code
299      !!-------------------------------------------------------------------
300      NAMELIST/namicetrp/ bound
301      !!-------------------------------------------------------------------
302
303      ! Read Namelist namicetrp
304      REWIND ( numnam_ice )
305      READ   ( numnam_ice  , namicetrp )
306      IF(lwp) THEN
307         WRITE(numout,*)
308         WRITE(numout,*) 'lim_trp_init_2 : Ice parameters for advection '
309         WRITE(numout,*) '~~~~~~~~~~~~~~'
310         WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound
311      ENDIF
312           
313   END SUBROUTINE lim_trp_init_2
314
315#else
316   !!----------------------------------------------------------------------
317   !!   Default option         Empty Module                No sea-ice model
318   !!----------------------------------------------------------------------
319CONTAINS
320   SUBROUTINE lim_trp_2        ! Empty routine
321   END SUBROUTINE lim_trp_2
322#endif
323
324   !!======================================================================
325END MODULE limtrp_2
Note: See TracBrowser for help on using the repository browser.