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

source: trunk/NEMO/LIM_SRC_3/limmsh.F90 @ 921

Last change on this file since 921 was 921, checked in by rblod, 16 years ago

Correct indentation and print for debug in LIM3, see ticket #134, step I

File size: 9.9 KB
Line 
1MODULE limmsh
2   !!======================================================================
3   !!                     ***  MODULE  limmsh  ***
4   !! LIM ice model :   definition of the ice mesh parameters
5   !!======================================================================
6#if defined key_lim3
7   !!----------------------------------------------------------------------
8   !!   'key_lim3'                                      LIM3 sea-ice model
9   !!----------------------------------------------------------------------
10   !!   lim_msh   : definition of the ice mesh
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE phycst
14   USE dom_oce
15   USE dom_ice
16   USE lbclnk
17   USE in_out_manager
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Accessibility
23   PUBLIC lim_msh      ! routine called by ice_ini.F90
24
25   !!----------------------------------------------------------------------
26   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
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
34      !!-------------------------------------------------------------------
35      !!                  ***  ROUTINE lim_msh  ***
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      !!---------------------------------------------------------------------
50      !! * Local variables
51      INTEGER :: ji, jj      ! dummy loop indices
52
53      REAL(wp), DIMENSION(jpi,jpj) ::  &
54         zd2d1 , zd1d2       ! Derivative of zh2 (resp. zh1) in the x direction
55      !                   ! (resp. y direction) (defined at the center)
56      REAL(wp) ::         &
57         zh1p  , zh2p   , &  ! Idem zh1, zh2 for the bottom left corner of the grid
58         zd2d1p, zd1d2p , &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid
59         zusden, zusden2     ! temporary scalars
60      !!---------------------------------------------------------------------
61
62      IF(lwp) THEN
63         WRITE(numout,*)
64         WRITE(numout,*) 'lim_msh : LIM sea-ice model, mesh initialization'
65         WRITE(numout,*) '~~~~~~~'
66      ENDIF
67
68      !----------------------------------------------------------                         
69      !    Initialization of local and some global (common) variables
70      !------------------------------------------------------------------
71
72      njeq   = INT( jpj / 2 )   !i bug mpp potentiel
73      njeqm1 = njeq - 1 
74
75      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor
76
77      IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN   ! local domain include both hemisphere
78         l_jeq = .TRUE.
79         njeq  = 1
80         DO WHILE ( njeq <= jpj .AND. fcor(1,njeq) < 0.e0 )
81            njeq = njeq + 1
82         END DO
83         IF(lwp ) WRITE(numout,*) '          the equator is inside the domain at about njeq = ', njeq
84      ELSEIF( fcor(1,1) < 0.e0 ) THEN
85         l_jeq = .FALSE.
86         njeq = jpj
87         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the southern hemisphere: njeq = ', njeq
88      ELSE
89         l_jeq = .FALSE.
90         njeq = 2
91         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the northern hemisphere: njeq = ', njeq
92      ENDIF
93
94      njeqm1 = njeq - 1
95
96
97      !   For each grid, definition of geometric tables
98      !------------------------------------------------------------------
99
100      !-------------------
101      ! Conventions :    !
102      !-------------------
103      !  indices 1 \ 2 <-> localisation in the 2 direction x \ y
104      !  3rd indice <-> localisation on the mesh :
105      !  0 = Centre ;  1 = corner W x(i-1/2) ; 2 = corner S y(j-1/2) ;
106      !  3 = corner SW x(i-1/2),y(j-1/2)
107      !-------------------
108      !!ibug ???
109      akappa(:,:,:,:) = 0.e0
110      wght(:,:,:,:) = 0.e0
111      alambd(:,:,:,:,:,:) = 0.e0
112      tmu(:,:) = 0.e0
113      tmv(:,:) = 0.0e0 ! CGrid EVP
114      !!i
115      ! metric coefficients for sea ice dynamic
116      !----------------------------------------
117      !                                                       ! akappa
118      DO jj = 2, jpj
119         zd1d2(:,jj) = e1v(:,jj) - e1v(:,jj-1)
120      END DO
121      CALL lbc_lnk( zd1d2, 'T', -1. )
122
123      DO ji = 2, jpi
124         zd2d1(ji,:) = e2u(ji,:) - e2u(ji-1,:)
125      END DO
126      CALL lbc_lnk( zd2d1, 'T', -1. )
127
128      akappa(:,:,1,1) =        1.0 / ( 2.0 * e1t(:,:) )
129      akappa(:,:,1,2) = zd1d2(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) )
130      akappa(:,:,2,1) = zd2d1(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) )
131      akappa(:,:,2,2) =        1.0 / ( 2.0 * e2t(:,:) )
132
133      !                                                      ! weights (wght)
134      DO jj = 2, jpj
135         DO ji = 2, jpi
136            zusden = 1. / (  ( e1t(ji,jj) + e1t(ji-1,jj  ) )   &
137               &           * ( e2t(ji,jj) + e2t(ji  ,jj-1) ) )
138            wght(ji,jj,1,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj  )
139            wght(ji,jj,1,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj-1)
140            wght(ji,jj,2,1) = zusden * e1t(ji-1,jj) * e2t(ji,jj  )
141            wght(ji,jj,2,2) = zusden * e1t(ji-1,jj) * e2t(ji,jj-1)
142         END DO
143      END DO
144      CALL lbc_lnk( wght(:,:,1,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V-point
145      CALL lbc_lnk( wght(:,:,1,2), 'I', 1. )      ! the value of wght at jpj is wrong
146      CALL lbc_lnk( wght(:,:,2,1), 'I', 1. )      ! but it is never used
147      CALL lbc_lnk( wght(:,:,2,2), 'I', 1. )
148
149      ! Coefficients for divergence of the stress tensor
150      !-------------------------------------------------
151
152      DO jj = 2, jpj
153         DO ji = 2, jpi
154            zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2)   &
155               &   + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   &
156               &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1)   &
157               &   + e1t(ji-1,jj-1) * wght(ji,jj,1,1)
158
159            zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2)   &
160               &   + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   &
161               &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1)   &
162               &   + e2t(ji-1,jj-1) * wght(ji,jj,1,1)
163
164            zusden = 1.0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 )
165            zusden2 = zusden * 2.0 
166
167            zd1d2p = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   )
168            zd2d1p = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   )
169
170            alambd(ji,jj,2,2,2,1) = zusden2 * e2t(ji  ,jj-1)
171            alambd(ji,jj,2,2,2,2) = zusden2 * e2t(ji  ,jj  )
172            alambd(ji,jj,2,2,1,1) = zusden2 * e2t(ji-1,jj-1)
173            alambd(ji,jj,2,2,1,2) = zusden2 * e2t(ji-1,jj  )
174
175            alambd(ji,jj,1,1,2,1) = zusden2 * e1t(ji  ,jj-1)
176            alambd(ji,jj,1,1,2,2) = zusden2 * e1t(ji  ,jj  )
177            alambd(ji,jj,1,1,1,1) = zusden2 * e1t(ji-1,jj-1)
178            alambd(ji,jj,1,1,1,2) = zusden2 * e1t(ji-1,jj  )
179
180            alambd(ji,jj,1,2,2,1) = zd1d2p
181            alambd(ji,jj,1,2,2,2) = zd1d2p
182            alambd(ji,jj,1,2,1,1) = zd1d2p
183            alambd(ji,jj,1,2,1,2) = zd1d2p
184
185            alambd(ji,jj,2,1,2,1) = zd2d1p
186            alambd(ji,jj,2,1,2,2) = zd2d1p
187            alambd(ji,jj,2,1,1,1) = zd2d1p
188            alambd(ji,jj,2,1,1,2) = zd2d1p
189         END DO
190      END DO
191
192      CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V point
193      CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. )      ! the value of wght at jpj is wrong
194      CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )      ! but it is never used
195      CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )      !
196
197      CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. )      ! CAUTION: idem
198      CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )      !
199      CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )      !
200      CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. )      !
201
202      CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. )      ! CAUTION: idem
203      CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. )      !
204      CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )      !
205      CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )      !
206
207      CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. )      ! CAUTION: idem
208      CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. )      !
209      CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )      !
210      CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )      !
211
212
213      ! Initialization of ice masks
214      !----------------------------
215
216      tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask
217
218      !      tmu(:,1) = 0.e0
219      !      tmu(1,:) = 0.e0
220      !      tmv(:,1) = 0.e0
221      !      tmv(1,:) = 0.e0
222
223      DO jj = 1, jpj - 1
224         DO ji = 1 , jpi - 1
225            tmu(ji,jj) =  tms(ji,jj) * tms(ji+1,jj)
226            tmv(ji,jj) =  tms(ji,jj) * tms(ji,jj+1)
227            tmf(ji,jj) =  tms(ji,jj) * tms(ji+1,jj) * tms(ji,jj+1) * &
228               tms(ji+1,jj+1)
229         END DO
230      END DO
231
232      !--lateral boundary conditions   
233      CALL lbc_lnk( tmu(:,:), 'U', 1. )
234      CALL lbc_lnk( tmv(:,:), 'V', 1. )
235      CALL lbc_lnk( tmf(:,:), 'F', 1. )
236
237      ! unmasked and masked area of T-grid cell
238      area(:,:) = e1t(:,:) * e2t(:,:)
239
240   END SUBROUTINE lim_msh
241
242#else
243   !!----------------------------------------------------------------------
244   !!   Default option            Dummy Module         NO LIM sea-ice model
245   !!----------------------------------------------------------------------
246CONTAINS
247   SUBROUTINE lim_msh           ! Dummy routine
248   END SUBROUTINE lim_msh
249#endif
250
251   !!======================================================================
252END MODULE limmsh
Note: See TracBrowser for help on using the repository browser.