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.F90 in trunk/NEMO/LIM_SRC – NEMO

source: trunk/NEMO/LIM_SRC/limtrp.F90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

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