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 @ 3

Last change on this file since 3 was 3, checked in by opalod, 21 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.6 KB
Line 
1MODULE limmsh
2#if defined key_ice_lim
3   !!======================================================================
4   !!                     ***  MODULE  limmsh  ***
5   !!                    definition of the ice mesh
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   lim_msh   : definition of the ice mesh
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE phycst
13   USE dom_oce
14   USE dom_ice
15   USE lbclnk
16
17   IMPLICIT NONE
18   PRIVATE
19
20   !! * Accessibility
21   PUBLIC lim_msh      ! routine called by ice_ini.F90
22
23   !!----------------------------------------------------------------------
24   !!   LIM 2.0 , UCL-LODYC-IPSL  (2003)
25   !!----------------------------------------------------------------------
26
27CONTAINS
28
29   SUBROUTINE lim_msh
30      !!-------------------------------------------------------------------
31      !!                  ***  ROUTINE lim_msh  ***
32      !!             
33      !! ** Purpose : Definition of the charact. of the numerical grid
34      !!       
35      !! ** Action  : - Initialisation of some variables
36      !!              - Definition of some constants linked with the grid
37      !!              - Definition of the metric coef. for the sea/ice
38      !!              - Initialization of the ice masks (tmsk, umsk)
39      !!
40      !! ** Refer.  : Deleersnijder et al. Ocean Modelling 100, 7-10
41      !!
42      !! ** History :
43      !!         original    : 01-04 (LIM)
44      !!         addition    : 02-08 (C. Ethe, G. Madec)
45      !!---------------------------------------------------------------------
46      !! * Local variables
47      INTEGER :: ji, jj     ! dummy loop indices
48
49      REAL(wp), DIMENSION(jpi,jpj) ::  &
50         zd2d1 , zd1d2      !  Derivative of zh2 (resp. zh1) in the x direction (defined at the center)
51
52      REAL(wp) ::         &
53         zh1p  , zh2p  ,  & !  Idem zh1, zh2 for the bottom left corner of the grid
54         zd2d1p, zd1d2p,  & !  Idem zd2d1, zd1d2 for the bottom left corner of the grid
55         zusden, zusden2, zaire4  !  tempory scalars
56      !!---------------------------------------------------------------------
57      !!  LIM 2.0, UCL-LODYC-IPSL (2002)
58      !!---------------------------------------------------------------------
59     
60      !----------------------------------------------------------                         
61      !    Initialization of local and some global (common) variables
62      !------------------------------------------------------------------
63     
64      jeq   = INT( jpj / 2 )   !i bug mpp potentiel
65      jeqm1 = jeq - 1 
66
67      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor
68 
69
70      !   For each grid, definition of geometric tables
71      !------------------------------------------------------------------
72     
73      !-------------------
74      ! Conventions :    |
75      !-------------------
76      !  indices 1 \ 2 <-> localisation in the 2 direction x \ y
77      !  3rd indice <-> localisation on the mesh :
78      !  0 = Centre ;  1 = corner W x(i-1/2) ; 2 = corner S y(j-1/2) ;
79      !  3 = corner SW x(i-1/2),y(j-1/2)
80      !-------------------
81     
82     
83      ! metric coefficients for sea ice dynamic
84      !----------------------------------------
85     
86      !                                                       ! akappa
87      DO jj = 2 , jpj
88         DO ji = 1 , jpi
89            zd1d2(ji,jj) = e1v(ji,jj) - e1v(ji,jj-1)
90         END DO
91      END DO
92      CALL lbc_lnk( zd1d2 , 'T', -1. )
93
94      DO jj = 1 , jpj
95         DO ji = 2 , jpi
96            zd2d1(ji,jj) = e2u(ji,jj) - e2u(ji-1,jj)
97         END DO
98      END DO
99      CALL lbc_lnk( zd2d1 , 'T', -1. )
100
101      DO jj = 1 , jpj
102         DO ji = 1 , jpi
103            zaire4            = 4.0 * e1t(ji,jj) * e2t(ji,jj)
104            akappa(ji,jj,1,1) =          1.0 / ( 2.0 * e1t(ji,jj) )
105            akappa(ji,jj,1,2) = zd1d2(ji,jj) / zaire4
106            akappa(ji,jj,2,1) = zd2d1(ji,jj) / zaire4
107            akappa(ji,jj,2,2) =          1.0 / ( 2.0 * e2t(ji,jj) )
108         END DO
109      END DO
110     
111      !                                                      ! weights (wght)
112      DO jj = 2 , jpj
113         DO ji = 2 , jpi
114            zusden = 1. / (  ( e1t(ji,jj) + e1t(ji-1,jj  ) )   &
115               &           * ( e2t(ji,jj) + e2t(ji  ,jj-1) ) )
116            wght(ji,jj,1,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj  )
117            wght(ji,jj,1,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj-1)
118            wght(ji,jj,2,1) = zusden * e1t(ji-1,jj) * e2t(ji,jj  )
119            wght(ji,jj,2,2) = zusden * e1t(ji-1,jj) * e2t(ji,jj-1)
120         END DO
121      END DO
122      CALL lbc_lnk( wght(:,:,1,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V-point
123      CALL lbc_lnk( wght(:,:,1,2), 'I', 1. )      ! the value of wght at jpj is wrong
124      CALL lbc_lnk( wght(:,:,2,1), 'I', 1. )      ! but it is never used
125      CALL lbc_lnk( wght(:,:,2,2), 'I', 1. )
126   
127      ! Coefficients for divergence of the stress tensor
128      !-------------------------------------------------
129
130      DO jj = 2,jpj
131         DO ji = 2 , jpi
132            zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2)   &
133               &   + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   &
134               &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1)   &
135               &   + e1t(ji-1,jj-1) * wght(ji,jj,1,1)
136
137            zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2)   &
138               &   + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   &
139               &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1)   &
140               &   + e2t(ji-1,jj-1) * wght(ji,jj,1,1)
141
142! better writen but change the last digit and thus solver in less than 100 timestep
143!           zh1p  = e1t(ji-1,jj  ) * wght(ji,jj,1,2) + e1t(ji,jj  ) * wght(ji,jj,2,2)   &
144!              &  + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1)
145
146!           zh2p  = e2t(ji-1,jj  ) * wght(ji,jj,1,2) + e2t(ji,jj  ) * wght(ji,jj,2,2)   &
147!              &  + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1)
148
149            zusden = 1.0 / ( zh1p * zh2p * 4.e0 )
150            zusden2 = zusden * 2.0 
151
152            zd1d2p = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   )
153            zd2d1p = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   )
154
155            alambd(ji,jj,2,2,2,1) = zusden2 * e2t(ji  ,jj-1)
156            alambd(ji,jj,2,2,2,2) = zusden2 * e2t(ji  ,jj  )
157            alambd(ji,jj,2,2,1,1) = zusden2 * e2t(ji-1,jj-1)
158            alambd(ji,jj,2,2,1,2) = zusden2 * e2t(ji-1,jj  )
159
160            alambd(ji,jj,1,1,2,1) = zusden2 * e1t(ji  ,jj-1)
161            alambd(ji,jj,1,1,2,2) = zusden2 * e1t(ji  ,jj  )
162            alambd(ji,jj,1,1,1,1) = zusden2 * e1t(ji-1,jj-1)
163            alambd(ji,jj,1,1,1,2) = zusden2 * e1t(ji-1,jj  )
164
165            alambd(ji,jj,1,2,2,1) = zd1d2p
166            alambd(ji,jj,1,2,2,2) = zd1d2p
167            alambd(ji,jj,1,2,1,1) = zd1d2p
168            alambd(ji,jj,1,2,1,2) = zd1d2p
169
170            alambd(ji,jj,2,1,2,1) = zd2d1p
171            alambd(ji,jj,2,1,2,2) = zd2d1p
172            alambd(ji,jj,2,1,1,1) = zd2d1p
173            alambd(ji,jj,2,1,1,2) = zd2d1p
174         END DO
175      END DO
176
177      CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V point
178      CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. )      ! the value of wght at jpj is wrong
179      CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )      ! but it is never used
180      CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )      !
181
182      CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. )      ! CAUTION: idem
183      CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )      !
184      CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )      !
185      CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. )      !
186
187      CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. )      ! CAUTION: idem
188      CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. )      !
189      CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )      !
190      CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )      !
191
192      CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. )      ! CAUTION: idem
193      CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. )      !
194      CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )      !
195      CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )      !
196           
197      ! Definition of scale dephts : bathymetry               
198      !-------------------------------------------
199!i bug  dz forced to 10 meters
200      dz = 10                  !!bug   potential bug if first level not equal to 10 m
201!!!      dz = gdept(1)
202
203     
204      ! Initialization of ice masks
205      !----------------------------
206     
207      tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask
208
209!i here we can use umask with a i and j shift of -1,-1
210      tmu(:,1) = 0.e0
211      tmu(1,:) = 0.e0
212      DO jj = 2, jpj               ! ice U.V-point: computed from ice T-point mask
213         DO ji = 2 , jpim1
214            tmu(ji,jj) =  tms(ji,jj) * tms(ji-1,jj) * tms(ji,jj-1) * tms(ji-1,jj-1)           
215         END DO
216      END DO
217     
218      !--lateral boundary conditions   
219      CALL lbc_lnk( tmu(:,:), 'I', 1. )
220     
221      ! unmasked and masked area of T-grid cell
222      area(:,:) = e1t(:,:) * e2t(:,:)
223      aire(:,:) = area(:,:) * tms(:,:)
224     
225   END SUBROUTINE lim_msh
226#else
227   !!==============================================================================
228   !!                       ***  MODULE limmsh   ***
229   !!                            No sea ice
230   !!==============================================================================
231CONTAINS
232   SUBROUTINE lim_msh           ! Empty routine
233   END SUBROUTINE lim_msh
234
235#endif
236   !!======================================================================
237END MODULE limmsh
Note: See TracBrowser for help on using the repository browser.