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

source: branches/dev_1784_EVP/NEMO/LIM_SRC_2/limtrp_2.F90 @ 2300

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

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