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

source: trunk/NEMO/LIM_SRC_2/limtrp_2.F90 @ 1922

Last change on this file since 1922 was 1922, checked in by rblod, 14 years ago

Correct alternate direction switch frenquence in LIM2 and LIM3, see ticket #625 and #605

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 16.1 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   !!----------------------------------------------------------------------
10#if defined key_lim2
11   !!----------------------------------------------------------------------
12   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_trp_2      : advection/diffusion process of sea ice
15   !!   lim_trp_init_2 : initialization and namelist read
16   !!----------------------------------------------------------------------
17   USE phycst          ! physical constant
18   USE sbc_oce         ! ocean surface boundary condition
19   USE dom_oce         ! ocean domain
20   USE in_out_manager  ! I/O manager
21   USE dom_ice_2       ! LIM-2 domain
22   USE ice_2           ! LIM-2 variables
23   USE limistate_2     ! LIM-2 initial state
24   USE limadv_2        ! LIM-2 advection
25   USE limhdf_2        ! LIM-2 horizontal diffusion
26   USE lbclnk          ! lateral boundary conditions -- 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)  ::           &  ! 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   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence (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      INTEGER, INTENT(in) ::   kt     ! number of iteration
66      !!
67      INTEGER  ::   ji, jj, jk   ! dummy loop indices
68      INTEGER  ::   initad       ! number of sub-timestep for the advection
69      REAL(wp) ::   zindb  , zindsn , zindic, zacrith   ! local scalars
70      REAL(wp) ::   zusvosn, zusvoic, zignm , zindhe    !   -      -
71      REAL(wp) ::   zvbord , zcfl   , zusnit            !   -      -
72      REAL(wp) ::   zrtt   , ztsn   , ztic1 , ztic2     !   -      -
73      REAL(wp), DIMENSION(jpi,jpj)  ::   zui_u , zvi_v , zsm             ! 2D workspace
74      REAL(wp), DIMENSION(jpi,jpj)  ::   zs0ice, zs0sn , zs0a            !  -      -
75      REAL(wp), DIMENSION(jpi,jpj)  ::   zs0c0 , zs0c1 , zs0c2 , zs0st   !  -      -
76      !---------------------------------------------------------------------
77
78      IF( kt == nit000  )   CALL lim_trp_init_2      ! Initialization (first time-step only)
79
80      zsm(:,:) = area(:,:)
81     
82      IF( ln_limdyn ) THEN
83         !-------------------------------------!
84         !   Advection of sea ice properties   !
85         !-------------------------------------!
86
87         ! ice velocities at ocean U- and V-points (zui_u,zvi_v)
88         ! ---------------------------------------
89         zvbord = 1.0 + ( 1.0 - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions       
90         DO jj = 1, jpjm1
91            DO ji = 1, jpim1   ! NO vector opt.
92               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 ) )
93               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 ) )
94            END DO
95         END DO
96         CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )         ! Lateral boundary conditions
97
98
99         ! CFL test for stability
100         ! ----------------------
101         zcfl  = 0.e0
102         zcfl  = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, :     ) ) * rdt_ice / e1u(1:jpim1, :     ) ) )
103         zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) )
104         !
105         IF(lk_mpp)   CALL mpp_max( zcfl )
106         !
107         IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ', zcfl
108
109         ! content of properties
110         ! ---------------------
111         zs0sn (:,:) =  hsnm(:,:) * area(:,:)                 ! Snow volume.
112         zs0ice(:,:) =  hicm(:,:) * area(:,:)                 ! Ice volume.
113         zs0a  (:,:) =  ( 1.0 - frld(:,:) )    * area  (:,:)  ! Surface covered by ice.
114         zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn (:,:)  ! Heat content of the snow layer.
115         zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)  ! Heat content of the first ice layer.
116         zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)  ! Heat content of the second ice layer.
117         zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a  (:,:)  ! Heat reservoir for brine pockets.
118         
119 
120         ! Advection (Prather scheme)
121         ! ---------
122         initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )   ! If ice drift field is too fast,         
123         zusnit = 1.0 / REAL( initad )                              ! split the ice time step in two
124         !
125         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN        !==  odd ice time step:  adv_x then adv_y  ==!
126            DO jk = 1, initad
127               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
128               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
129               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
130               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
131               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
132               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
133               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
134               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
135               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
136               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
137               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
138               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
139               CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
140               CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
141            END DO
142         ELSE                                                 !==  even ice time step:  adv_x then adv_y  ==!
143            DO jk = 1, initad
144               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
145               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
146               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
147               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn  )
148               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
149               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0a  , sxa  , sxxa  , sya  , syya  , sxya   )
150               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
151               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0  )
152               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
153               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1  )
154               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
155               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2  )
156               CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
157               CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst  )
158            END DO
159         ENDIF
160                       
161         ! recover the properties from their contents
162         ! ------------------------------------------
163!!gm Define in limmsh one for all area = 1 /area  (CPU time saved !)
164         zs0ice(:,:) = zs0ice(:,:) / area(:,:)
165         zs0sn (:,:) = zs0sn (:,:) / area(:,:)
166         zs0a  (:,:) = zs0a  (:,:) / area(:,:)
167         zs0c0 (:,:) = zs0c0 (:,:) / area(:,:)
168         zs0c1 (:,:) = zs0c1 (:,:) / area(:,:)
169         zs0c2 (:,:) = zs0c2 (:,:) / area(:,:)
170         zs0st (:,:) = zs0st (:,:) / area(:,:)
171
172
173         !-------------------------------------!
174         !   Diffusion of sea ice properties   !
175         !-------------------------------------!
176
177         ! Masked eddy diffusivity coefficient at ocean U- and V-points
178         ! ------------------------------------------------------------
179         DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
180            DO ji = 1 , fs_jpim1   ! vector opt.
181               pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji  ,jj) ) ) )   &
182                  &        * ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj) ) ) ) * ahiu(ji,jj)
183               pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji,jj  ) ) ) )   &
184                  &        * ( 1.0 - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1) ) ) ) * ahiv(ji,jj)
185            END DO
186         END DO
187!!gm more readable coding: (and avoid an error in F90 with sign of zero)
188!        DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
189!           DO ji = 1 , fs_jpim1   ! vector opt.
190!              IF( MIN( zs0a(ji,jj) , zs0a(ji+1,jj) ) == 0.e0 )   pahu(ji,jj) = 0.e0
191!              IF( MIN( zs0a(ji,jj) , zs0a(ji,jj+1) ) == 0.e0 )   pahv(ji,jj) = 0.e0
192!           END DO
193!        END DO
194!!gm end
195
196         ! diffusion
197         ! ---------
198         CALL lim_hdf_2( zs0ice )
199         CALL lim_hdf_2( zs0sn  )
200         CALL lim_hdf_2( zs0a   )
201         CALL lim_hdf_2( zs0c0  )
202         CALL lim_hdf_2( zs0c1  )
203         CALL lim_hdf_2( zs0c2  )
204         CALL lim_hdf_2( zs0st  )
205
206!!gm see comment this can be skipped
207         zs0ice(:,:) = MAX( rzero, zs0ice(:,:) * area(:,:) )    !!bug:  useless
208         zs0sn (:,:) = MAX( rzero, zs0sn (:,:) * area(:,:) )    !!bug:  cf /area  just below
209         zs0a  (:,:) = MAX( rzero, zs0a  (:,:) * area(:,:) )    !! caution: the suppression of the 2 changes
210         zs0c0 (:,:) = MAX( rzero, zs0c0 (:,:) * area(:,:) )    !! the last digit of the results
211         zs0c1 (:,:) = MAX( rzero, zs0c1 (:,:) * area(:,:) )
212         zs0c2 (:,:) = MAX( rzero, zs0c2 (:,:) * area(:,:) )
213         zs0st (:,:) = MAX( rzero, zs0st (:,:) * area(:,:) )
214
215
216         !-------------------------------------------------------------------!
217         !   Updating and limitation of sea ice properties after transport   !
218         !-------------------------------------------------------------------!
219         DO jj = 1, jpj
220            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH
221            DO ji = 1, jpi
222               !
223               ! Recover mean values over the grid squares.
224               zs0sn (ji,jj) = MAX( rzero, zs0sn (ji,jj)/area(ji,jj) )
225               zs0ice(ji,jj) = MAX( rzero, zs0ice(ji,jj)/area(ji,jj) )
226               zs0a  (ji,jj) = MAX( rzero, zs0a  (ji,jj)/area(ji,jj) )
227               zs0c0 (ji,jj) = MAX( rzero, zs0c0 (ji,jj)/area(ji,jj) )
228               zs0c1 (ji,jj) = MAX( rzero, zs0c1 (ji,jj)/area(ji,jj) )
229               zs0c2 (ji,jj) = MAX( rzero, zs0c2 (ji,jj)/area(ji,jj) )
230               zs0st (ji,jj) = MAX( rzero, zs0st (ji,jj)/area(ji,jj) )
231
232               ! Recover in situ values.
233               zindb         = MAX( rzero, SIGN( rone, zs0a(ji,jj) - epsi06 ) )
234               zacrith       = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) )
235               zs0a (ji,jj)  = zindb * MIN( zs0a(ji,jj), zacrith )
236               hsnif(ji,jj)  = zindb * ( zs0sn(ji,jj) /MAX( zs0a(ji,jj), epsi16 ) )
237               hicif(ji,jj)  = zindb * ( zs0ice(ji,jj)/MAX( zs0a(ji,jj), epsi16 ) )
238               zindsn        = MAX( rzero, SIGN( rone, hsnif(ji,jj) - epsi06 ) )
239               zindic        = MAX( rzero, SIGN( rone, hicif(ji,jj) - epsi03 ) )
240               zindb         = MAX( zindsn, zindic )
241               zs0a (ji,jj)  = zindb * zs0a(ji,jj)
242               frld (ji,jj)  = 1.0 - zs0a(ji,jj)
243               hsnif(ji,jj)  = zindsn * hsnif(ji,jj)
244               hicif(ji,jj)  = zindic * hicif(ji,jj)
245               zusvosn       = 1.0/MAX( hsnif(ji,jj) * zs0a(ji,jj), epsi16 )
246               zusvoic       = 1.0/MAX( hicif(ji,jj) * zs0a(ji,jj), epsi16 )
247               zignm         = MAX( rzero,  SIGN( rone, hsndif - hsnif(ji,jj) ) )
248               zrtt          = 173.15 * rone 
249               ztsn          =          zignm   * tbif(ji,jj,1)  &
250                              + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) ) 
251               ztic1          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c1(ji,jj) ) , tfu(ji,jj) )
252               ztic2          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c2(ji,jj) ) , tfu(ji,jj) )
253 
254               tbif(ji,jj,1) = zindsn * ztsn  + ( 1.0 - zindsn ) * tfu(ji,jj)               
255               tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj)
256               tbif(ji,jj,3) = zindic * ztic2 + ( 1.0 - zindic ) * tfu(ji,jj)
257               qstoif(ji,jj) = zindb  * xlic * zs0st(ji,jj) /  MAX( zs0a(ji,jj), epsi16 )
258            END DO
259         END DO
260         !
261      ENDIF
262      !
263   END SUBROUTINE lim_trp_2
264
265
266   SUBROUTINE lim_trp_init_2
267      !!-------------------------------------------------------------------
268      !!                  ***  ROUTINE lim_trp_init_2  ***
269      !!
270      !! ** Purpose :   initialization of ice advection parameters
271      !!
272      !! ** Method  :   Read the namicetrp namelist and check the parameter
273      !!              values called at the first timestep (nit000)
274      !!
275      !! ** input   :   Namelist namicetrp
276      !!-------------------------------------------------------------------
277      NAMELIST/namicetrp/ bound
278      !!-------------------------------------------------------------------
279      !
280      REWIND ( numnam_ice )      ! Read Namelist namicetrp
281      READ   ( numnam_ice  , namicetrp )
282      IF(lwp) THEN
283         WRITE(numout,*)
284         WRITE(numout,*) 'lim_trp_init_2 : Ice parameters for advection '
285         WRITE(numout,*) '~~~~~~~~~~~~~~'
286         WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound
287      ENDIF
288      !
289   END SUBROUTINE lim_trp_init_2
290
291#else
292   !!----------------------------------------------------------------------
293   !!   Default option         Empty Module                No sea-ice model
294   !!----------------------------------------------------------------------
295CONTAINS
296   SUBROUTINE lim_trp_2        ! Empty routine
297   END SUBROUTINE lim_trp_2
298#endif
299
300   !!======================================================================
301END MODULE limtrp_2
Note: See TracBrowser for help on using the repository browser.