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

source: branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limtrp_2.F90 @ 1855

Last change on this file since 1855 was 1855, checked in by gm, 14 years ago

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