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 – NEMO

source: trunk/NEMO/LIM_SRC/limmsh.F90 @ 192

Last change on this file since 192 was 192, checked in by opalod, 19 years ago

CT : BUGFIX133 : correction of the njeq initialisation to avoid compilation error when making a bound array check
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.4 KB
Line 
1MODULE limmsh
2   !!======================================================================
3   !!                     ***  MODULE  limmsh  ***
4   !! LIM ice model :   definition of the ice mesh parameters
5   !!======================================================================
6#if defined key_ice_lim
7   !!----------------------------------------------------------------------
8   !!   'key_ice_lim'                                     LIM 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 2.0 , UCL-LODYC-IPSL  (2003)
27   !!----------------------------------------------------------------------
28
29CONTAINS
30
31   SUBROUTINE lim_msh
32      !!-------------------------------------------------------------------
33      !!                  ***  ROUTINE lim_msh  ***
34      !!             
35      !! ** Purpose : Definition of the charact. of the numerical grid
36      !!       
37      !! ** Action  : - Initialisation of some variables
38      !!              - Definition of some constants linked with the grid
39      !!              - Definition of the metric coef. for the sea/ice
40      !!              - Initialization of the ice masks (tmsk, umsk)
41      !!
42      !! ** Refer.  : Deleersnijder et al. Ocean Modelling 100, 7-10
43      !!
44      !! ** History :
45      !!         original    : 01-04 (LIM)
46      !!         addition    : 02-08 (C. Ethe, G. Madec)
47      !!---------------------------------------------------------------------
48      !! * Local variables
49      INTEGER :: ji, jj      ! dummy loop indices
50
51      REAL(wp), DIMENSION(jpi,jpj) ::  &
52         zd2d1 , zd1d2       ! Derivative of zh2 (resp. zh1) in the x direction
53         !                   ! (resp. y direction) (defined at the center)
54      REAL(wp) ::         &
55         zh1p  , zh2p   , &  ! Idem zh1, zh2 for the bottom left corner of the grid
56         zd2d1p, zd1d2p , &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid
57         zusden, zusden2     ! temporary scalars
58      !!---------------------------------------------------------------------
59
60      IF(lwp) THEN
61         WRITE(numout,*)
62         WRITE(numout,*) 'lim_msh : LIM sea-ice model, mesh initialization'
63         WRITE(numout,*) '~~~~~~~'
64      ENDIF
65     
66      !----------------------------------------------------------                         
67      !    Initialization of local and some global (common) variables
68      !------------------------------------------------------------------
69     
70      njeq   = INT( jpj / 2 )   !i bug mpp potentiel
71      njeqm1 = njeq - 1 
72
73      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor
74 
75!i    DO jj = 1, jpj
76!i       zmsk(jj) = SUM( tmask(:,jj,:) )   ! = 0          if land  everywhere on a j-line
77!!ii     write(numout,*) jj, zind(jj)
78!i    END DO
79
80      IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN   ! local domain include both hemisphere
81         l_jeq = .TRUE.
82         njeq  = 1
83         DO WHILE ( njeq <= jpj .AND. fcor(1,njeq) < 0.e0 )
84            njeq = njeq + 1
85         END DO
86         IF(lwp ) WRITE(numout,*) '          the equator is inside the domain at about njeq = ', njeq
87      ELSEIF( fcor(1,1) < 0.e0 ) THEN
88         l_jeq = .FALSE.
89         njeq = jpj
90         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the southern hemisphere: njeq = ', njeq
91      ELSE
92         l_jeq = .FALSE.
93         njeq = 2
94         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the northern hemisphere: njeq = ', njeq
95      ENDIF
96
97      njeqm1 = njeq - 1
98
99
100      !   For each grid, definition of geometric tables
101      !------------------------------------------------------------------
102     
103      !-------------------
104      ! Conventions :    !
105      !-------------------
106      !  indices 1 \ 2 <-> localisation in the 2 direction x \ y
107      !  3rd indice <-> localisation on the mesh :
108      !  0 = Centre ;  1 = corner W x(i-1/2) ; 2 = corner S y(j-1/2) ;
109      !  3 = corner SW x(i-1/2),y(j-1/2)
110      !-------------------
111!!ibug ???
112      akappa(:,:,:,:) = 0.e0
113      wght(:,:,:,:) = 0.e0
114      alambd(:,:,:,:,:,:) = 0.e0
115      tmu(:,:) = 0.e0
116!!i
117     
118     
119      ! metric coefficients for sea ice dynamic
120      !----------------------------------------
121      !                                                       ! akappa
122      DO jj = 2, jpj
123         zd1d2(:,jj) = e1v(:,jj) - e1v(:,jj-1)
124      END DO
125      CALL lbc_lnk( zd1d2, 'T', -1. )
126
127      DO ji = 2, jpi
128         zd2d1(ji,:) = e2u(ji,:) - e2u(ji-1,:)
129      END DO
130      CALL lbc_lnk( zd2d1, 'T', -1. )
131
132      akappa(:,:,1,1) =        1.0 / ( 2.0 * e1t(:,:) )
133      akappa(:,:,1,2) = zd1d2(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) )
134      akappa(:,:,2,1) = zd2d1(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) )
135      akappa(:,:,2,2) =        1.0 / ( 2.0 * e2t(:,:) )
136     
137      !                                                      ! weights (wght)
138      DO jj = 2, jpj
139         DO ji = 2, jpi
140            zusden = 1. / (  ( e1t(ji,jj) + e1t(ji-1,jj  ) )   &
141               &           * ( e2t(ji,jj) + e2t(ji  ,jj-1) ) )
142            wght(ji,jj,1,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj  )
143            wght(ji,jj,1,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj-1)
144            wght(ji,jj,2,1) = zusden * e1t(ji-1,jj) * e2t(ji,jj  )
145            wght(ji,jj,2,2) = zusden * e1t(ji-1,jj) * e2t(ji,jj-1)
146         END DO
147      END DO
148      CALL lbc_lnk( wght(:,:,1,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V-point
149      CALL lbc_lnk( wght(:,:,1,2), 'I', 1. )      ! the value of wght at jpj is wrong
150      CALL lbc_lnk( wght(:,:,2,1), 'I', 1. )      ! but it is never used
151      CALL lbc_lnk( wght(:,:,2,2), 'I', 1. )
152   
153      ! Coefficients for divergence of the stress tensor
154      !-------------------------------------------------
155
156      DO jj = 2, jpj
157         DO ji = 2, jpi
158            zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2)   &
159               &   + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   &
160               &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1)   &
161               &   + e1t(ji-1,jj-1) * wght(ji,jj,1,1)
162
163            zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2)   &
164               &   + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   &
165               &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1)   &
166               &   + e2t(ji-1,jj-1) * wght(ji,jj,1,1)
167
168! better written but change the last digit and thus solver in less than 100 timestep
169!           zh1p  = e1t(ji-1,jj  ) * wght(ji,jj,1,2) + e1t(ji,jj  ) * wght(ji,jj,2,2)   &
170!              &  + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1)
171
172!           zh2p  = e2t(ji-1,jj  ) * wght(ji,jj,1,2) + e2t(ji,jj  ) * wght(ji,jj,2,2)   &
173!              &  + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1)
174
175!!ibug =0   zusden = 1.0 / ( zh1p * zh2p * 4.e0 )
176            zusden = 1.0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 )
177            zusden2 = zusden * 2.0 
178
179            zd1d2p = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   )
180            zd2d1p = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   )
181
182            alambd(ji,jj,2,2,2,1) = zusden2 * e2t(ji  ,jj-1)
183            alambd(ji,jj,2,2,2,2) = zusden2 * e2t(ji  ,jj  )
184            alambd(ji,jj,2,2,1,1) = zusden2 * e2t(ji-1,jj-1)
185            alambd(ji,jj,2,2,1,2) = zusden2 * e2t(ji-1,jj  )
186
187            alambd(ji,jj,1,1,2,1) = zusden2 * e1t(ji  ,jj-1)
188            alambd(ji,jj,1,1,2,2) = zusden2 * e1t(ji  ,jj  )
189            alambd(ji,jj,1,1,1,1) = zusden2 * e1t(ji-1,jj-1)
190            alambd(ji,jj,1,1,1,2) = zusden2 * e1t(ji-1,jj  )
191
192            alambd(ji,jj,1,2,2,1) = zd1d2p
193            alambd(ji,jj,1,2,2,2) = zd1d2p
194            alambd(ji,jj,1,2,1,1) = zd1d2p
195            alambd(ji,jj,1,2,1,2) = zd1d2p
196
197            alambd(ji,jj,2,1,2,1) = zd2d1p
198            alambd(ji,jj,2,1,2,2) = zd2d1p
199            alambd(ji,jj,2,1,1,1) = zd2d1p
200            alambd(ji,jj,2,1,1,2) = zd2d1p
201         END DO
202      END DO
203
204      CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V point
205      CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. )      ! the value of wght at jpj is wrong
206      CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )      ! but it is never used
207      CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )      !
208
209      CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. )      ! CAUTION: idem
210      CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )      !
211      CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )      !
212      CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. )      !
213
214      CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. )      ! CAUTION: idem
215      CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. )      !
216      CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )      !
217      CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )      !
218
219      CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. )      ! CAUTION: idem
220      CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. )      !
221      CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )      !
222      CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )      !
223           
224
225      ! Initialization of ice masks
226      !----------------------------
227     
228      tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask
229
230!i here we can use umask with a i and j shift of -1,-1
231      tmu(:,1) = 0.e0
232      tmu(1,:) = 0.e0
233      DO jj = 2, jpj               ! ice U.V-point: computed from ice T-point mask
234         DO ji = 2, jpim1
235            tmu(ji,jj) =  tms(ji,jj) * tms(ji-1,jj) * tms(ji,jj-1) * tms(ji-1,jj-1)           
236         END DO
237      END DO
238     
239      !--lateral boundary conditions   
240      CALL lbc_lnk( tmu(:,:), 'I', 1. )
241     
242      ! unmasked and masked area of T-grid cell
243      area(:,:) = e1t(:,:) * e2t(:,:)
244     
245   END SUBROUTINE lim_msh
246
247#else
248   !!----------------------------------------------------------------------
249   !!   Default option            Dummy Module         NO LIM sea-ice model
250   !!----------------------------------------------------------------------
251CONTAINS
252   SUBROUTINE lim_msh           ! Dummy routine
253   END SUBROUTINE lim_msh
254#endif
255
256   !!======================================================================
257END MODULE limmsh
Note: See TracBrowser for help on using the repository browser.