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.
limmsh_2.F90 in branches/devmercator2010/NEMO/LIM_SRC_2 – NEMO

source: branches/devmercator2010/NEMO/LIM_SRC_2/limmsh_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: 12.3 KB
Line 
1MODULE limmsh_2
2   !!======================================================================
3   !!                     ***  MODULE  limmsh_2  ***
4   !! LIM 2.0 ice model :   definition of the ice mesh parameters
5   !!======================================================================
6#if defined key_lim2
7   !!----------------------------------------------------------------------
8   !!   'key_lim2'                                     LIM 2.0sea-ice model
9   !!----------------------------------------------------------------------
10   !!   lim_msh_2   : definition of the ice mesh
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE phycst
14   USE dom_oce
15   USE dom_ice_2
16   USE lbclnk
17   USE in_out_manager
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Accessibility
23   PUBLIC lim_msh_2      ! routine called by ice_ini_2.F90
24
25   !!----------------------------------------------------------------------
26   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)
27   !! $Id$
28   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
29   !!----------------------------------------------------------------------
30
31CONTAINS
32
33   SUBROUTINE lim_msh_2
34      !!-------------------------------------------------------------------
35      !!                  ***  ROUTINE lim_msh_2  ***
36      !!             
37      !! ** Purpose : Definition of the charact. of the numerical grid
38      !!       
39      !! ** Action  : - Initialisation of some variables
40      !!              - Definition of some constants linked with the grid
41      !!              - Definition of the metric coef. for the sea/ice
42      !!              - Initialization of the ice masks (tmsk, umsk)
43      !!
44      !! ** Refer.  : Deleersnijder et al. Ocean Modelling 100, 7-10
45      !!
46      !! ** History :
47      !!         original    : 01-04 (LIM)
48      !!         addition    : 02-08 (C. Ethe, G. Madec)
49      !!         additions   : 2009-05 (addition of the lim2_evp case, G. Garric)
50      !!---------------------------------------------------------------------
51      !! * Local variables
52      INTEGER :: ji, jj      ! dummy loop indices
53
54      REAL(wp) ::         &
55         zusden              ! temporary scalars
56#if defined key_lim2_vp
57      REAL(wp), DIMENSION(jpi,jpj) ::  &
58         zd2d1 , zd1d2       ! Derivative of zh2 (resp. zh1) in the x direction
59         !                   ! (resp. y direction) (defined at the center)
60      REAL(wp) ::         &
61         zh1p  , zh2p   , &  ! Idem zh1, zh2 for the bottom left corner of the grid
62         zd2d1p, zd1d2p , &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid
63         zusden2             ! temporary scalars
64#endif
65      !!---------------------------------------------------------------------
66
67      IF(lwp) THEN
68         WRITE(numout,*)
69         WRITE(numout,*) 'lim_msh_2 : LIM 2.0 sea-ice model, mesh initialization'
70         WRITE(numout,*) '~~~~~~~~~'
71      ENDIF
72     
73      !----------------------------------------------------------                         
74      !    Initialization of local and some global (common) variables
75      !------------------------------------------------------------------
76     
77      njeq   = INT( jpj / 2 )   !i bug mpp potentiel
78      njeqm1 = njeq - 1 
79
80      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor
81 
82!i    DO jj = 1, jpj
83!i       zmsk(jj) = SUM( tmask(:,jj,:) )   ! = 0          if land  everywhere on a j-line
84!!ii     write(numout,*) jj, zind(jj)
85!i    END DO
86
87      IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN   ! local domain include both hemisphere
88         l_jeq = .TRUE.
89         njeq  = 1
90         DO WHILE ( njeq <= jpj .AND. fcor(1,njeq) < 0.e0 )
91            njeq = njeq + 1
92         END DO
93         IF(lwp ) WRITE(numout,*) '          the equator is inside the domain at about njeq = ', njeq
94      ELSEIF( fcor(1,1) < 0.e0 ) THEN
95         l_jeq = .FALSE.
96         njeq = jpj
97         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the southern hemisphere: njeq = ', njeq
98      ELSE
99         l_jeq = .FALSE.
100         njeq = 2
101         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the northern hemisphere: njeq = ', njeq
102      ENDIF
103
104      njeqm1 = njeq - 1
105
106
107      !   For each grid, definition of geometric tables
108      !------------------------------------------------------------------
109     
110      !-------------------
111      ! Conventions :    !
112      !-------------------
113      !  indices 1 \ 2 <-> localisation in the 2 direction x \ y
114      !  3rd indice <-> localisation on the mesh :
115      !  0 = Centre ;  1 = corner W x(i-1/2) ; 2 = corner S y(j-1/2) ;
116      !  3 = corner SW x(i-1/2),y(j-1/2)
117      !-------------------
118!!ibug ???
119      wght(:,:,:,:) = 0.e0
120      tmu(:,:)      = 0.e0
121#if defined key_lim2_vp 
122      akappa(:,:,:,:)     = 0.e0
123      alambd(:,:,:,:,:,:) = 0.e0
124#else
125      tmv(:,:) = 0.e0
126      tmf(:,:) = 0.e0
127#endif
128!!i
129     
130#if defined key_lim2_vp
131      ! metric coefficients for sea ice dynamic
132      !----------------------------------------
133      !                                                       ! akappa
134      DO jj = 2, jpj
135         zd1d2(:,jj) = e1v(:,jj) - e1v(:,jj-1)
136      END DO
137      CALL lbc_lnk( zd1d2, 'T', -1. )
138
139      DO ji = 2, jpi
140         zd2d1(ji,:) = e2u(ji,:) - e2u(ji-1,:)
141      END DO
142      CALL lbc_lnk( zd2d1, 'T', -1. )
143
144      akappa(:,:,1,1) =        1.0 / ( 2.0 * e1t(:,:) )
145      akappa(:,:,1,2) = zd1d2(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) )
146      akappa(:,:,2,1) = zd2d1(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) )
147      akappa(:,:,2,2) =        1.0 / ( 2.0 * e2t(:,:) )
148     
149      !                                                      ! weights (wght)
150      DO jj = 2, jpj
151         DO ji = 2, jpi
152            zusden = 1. / (  ( e1t(ji,jj) + e1t(ji-1,jj  ) )   &
153               &           * ( e2t(ji,jj) + e2t(ji  ,jj-1) ) )
154            wght(ji,jj,1,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj  )
155            wght(ji,jj,1,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj-1)
156            wght(ji,jj,2,1) = zusden * e1t(ji-1,jj) * e2t(ji,jj  )
157            wght(ji,jj,2,2) = zusden * e1t(ji-1,jj) * e2t(ji,jj-1)
158         END DO
159      END DO
160      CALL lbc_lnk( wght(:,:,1,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V-point
161      CALL lbc_lnk( wght(:,:,1,2), 'I', 1. )      ! the value of wght at jpj is wrong
162      CALL lbc_lnk( wght(:,:,2,1), 'I', 1. )      ! but it is never used
163      CALL lbc_lnk( wght(:,:,2,2), 'I', 1. )
164#else
165      !                                                      ! weights (wght)
166      DO jj = 2, jpj-1
167         DO ji = 2, jpi-1
168            zusden = 1. / (  ( e1t(ji+1,jj) + e1t(ji,jj  ) )   &
169               &           * ( e2t(ji,jj+1) + e2t(ji  ,jj) ) ) 
170            wght(ji,jj,1,1) = zusden * e1t(ji+1,jj) * e2t(ji,jj+1)
171            wght(ji,jj,1,2) = zusden * e1t(ji+1,jj) * e2t(ji,jj  )
172            wght(ji,jj,2,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj+1)
173            wght(ji,jj,2,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj  )
174         END DO
175      END DO
176
177      !With EVP, the weights are calculated on 'F' points
178      CALL lbc_lnk( wght(:,:,1,1), 'F', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V-point
179      CALL lbc_lnk( wght(:,:,1,2), 'F', 1. )      ! the value of wght at jpj is wrong
180      CALL lbc_lnk( wght(:,:,2,1), 'F', 1. )      ! but it is never used
181      CALL lbc_lnk( wght(:,:,2,2), 'F', 1. )
182
183#endif
184   
185      ! Coefficients for divergence of the stress tensor
186      !-------------------------------------------------
187
188#if defined key_lim2_vp
189      DO jj = 2, jpj
190         DO ji = 2, jpi   ! NO vector opt.
191            zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2)   &
192               &   + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   &
193               &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1)   &
194               &   + e1t(ji-1,jj-1) * wght(ji,jj,1,1)
195
196            zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2)   &
197               &   + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   &
198               &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1)   &
199               &   + e2t(ji-1,jj-1) * wght(ji,jj,1,1)
200
201! better written but change the last digit and thus solver in less than 100 timestep
202!           zh1p  = e1t(ji-1,jj  ) * wght(ji,jj,1,2) + e1t(ji,jj  ) * wght(ji,jj,2,2)   &
203!              &  + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1)
204
205!           zh2p  = e2t(ji-1,jj  ) * wght(ji,jj,1,2) + e2t(ji,jj  ) * wght(ji,jj,2,2)   &
206!              &  + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1)
207
208!!ibug =0   zusden = 1.0 / ( zh1p * zh2p * 4.e0 )
209            zusden = 1.0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 )
210            zusden2 = zusden * 2.0 
211
212            zd1d2p = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   )
213            zd2d1p = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   )
214
215            alambd(ji,jj,2,2,2,1) = zusden2 * e2t(ji  ,jj-1)
216            alambd(ji,jj,2,2,2,2) = zusden2 * e2t(ji  ,jj  )
217            alambd(ji,jj,2,2,1,1) = zusden2 * e2t(ji-1,jj-1)
218            alambd(ji,jj,2,2,1,2) = zusden2 * e2t(ji-1,jj  )
219
220            alambd(ji,jj,1,1,2,1) = zusden2 * e1t(ji  ,jj-1)
221            alambd(ji,jj,1,1,2,2) = zusden2 * e1t(ji  ,jj  )
222            alambd(ji,jj,1,1,1,1) = zusden2 * e1t(ji-1,jj-1)
223            alambd(ji,jj,1,1,1,2) = zusden2 * e1t(ji-1,jj  )
224
225            alambd(ji,jj,1,2,2,1) = zd1d2p
226            alambd(ji,jj,1,2,2,2) = zd1d2p
227            alambd(ji,jj,1,2,1,1) = zd1d2p
228            alambd(ji,jj,1,2,1,2) = zd1d2p
229
230            alambd(ji,jj,2,1,2,1) = zd2d1p
231            alambd(ji,jj,2,1,2,2) = zd2d1p
232            alambd(ji,jj,2,1,1,1) = zd2d1p
233            alambd(ji,jj,2,1,1,2) = zd2d1p
234         END DO
235      END DO
236
237      CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V point
238      CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. )      ! the value of wght at jpj is wrong
239      CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )      ! but it is never used
240      CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )      !
241
242      CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. )      ! CAUTION: idem
243      CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )      !
244      CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )      !
245      CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. )      !
246
247      CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. )      ! CAUTION: idem
248      CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. )      !
249      CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )      !
250      CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )      !
251
252      CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. )      ! CAUTION: idem
253      CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. )      !
254      CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )      !
255      CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )      !
256#endif
257           
258
259      ! Initialization of ice masks
260      !----------------------------
261     
262      tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask
263
264!i here we can use umask with a i and j shift of -1,-1
265      tmu(:,1) = 0.e0
266      tmu(1,:) = 0.e0
267
268#if defined key_lim2_vp
269      DO jj = 2, jpj               ! ice U.V-point: computed from ice T-point mask
270         DO ji = 2, jpim1   ! NO vector opt.
271            tmu(ji,jj) =  tms(ji,jj) * tms(ji-1,jj) * tms(ji,jj-1) * tms(ji-1,jj-1)           
272         END DO
273      END DO
274     
275      !--lateral boundary conditions   
276      CALL lbc_lnk( tmu(:,:), 'I', 1. )
277#else
278      tmv(:,1) = 0.e0 !SB
279      tmv(1,:) = 0.e0 !SB
280      tmf(1,:) = 0.e0
281      tmf(:,1) = 0.e0
282      DO jj = 1, jpj - 1
283         DO ji = 1 , jpi - 1
284            tmu(ji,jj) =  tms(ji,jj) * tms(ji+1,jj)
285            tmv(ji,jj) =  tms(ji,jj) * tms(ji,jj+1)
286            tmf(ji,jj) =  tms(ji,jj) * tms(ji+1,jj) * tms(ji,jj+1) * &
287               tms(ji+1,jj+1)
288         END DO
289      END DO
290
291      !--lateral boundary conditions
292      CALL lbc_lnk( tmu(:,:), 'U', 1. )
293      CALL lbc_lnk( tmv(:,:), 'V', 1. )
294      CALL lbc_lnk( tmf(:,:), 'F', 1. )
295#endif
296     
297      ! unmasked and masked area of T-grid cell
298      area(:,:) = e1t(:,:) * e2t(:,:)
299     
300   END SUBROUTINE lim_msh_2
301
302#else
303   !!----------------------------------------------------------------------
304   !!   Default option            Dummy Module         NO LIM sea-ice model
305   !!----------------------------------------------------------------------
306CONTAINS
307   SUBROUTINE lim_msh_2           ! Dummy routine
308   END SUBROUTINE lim_msh_2
309#endif
310
311   !!======================================================================
312END MODULE limmsh_2
Note: See TracBrowser for help on using the repository browser.