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

source: branches/dev_1784_EVP/NEMO/LIM_SRC_2/limmsh_2.F90 @ 2293

Last change on this file since 2293 was 2095, checked in by cbricaud, 14 years ago

add correction from reviewer
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 11.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   !! History :  LIM  ! 2001-04  (Louvain-la-Neuve) Original code
7   !!            1.0  ! 2002-08  (C. Ethe, G. Madec)
8   !!            3.3  ! 2009-05  (G. Garric, C. Bricaud) addition of the lim2_evp case
9   !!----------------------------------------------------------------------
10#if defined   key_lim2
11   !!----------------------------------------------------------------------
12   !!   'key_lim2'                                     LIM 2.0sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_msh_2   : definition of the ice mesh
15   !!----------------------------------------------------------------------
16   USE phycst           ! physical constants
17   USE dom_oce          ! ocean domain
18   USE dom_ice_2        ! LIM2: ice domain
19   USE lbclnk           ! lateral boundary condition - MPP exchanges
20   USE in_out_manager   ! I/O manager
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   lim_msh_2   ! routine called by ice_ini_2.F90
26
27   !!----------------------------------------------------------------------
28   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010)
29   !! $Id$
30   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
31   !!----------------------------------------------------------------------
32CONTAINS
33
34   SUBROUTINE lim_msh_2
35      !!-------------------------------------------------------------------
36      !!                  ***  ROUTINE lim_msh_2  ***
37      !!             
38      !! ** Purpose : Definition of the charact. of the numerical grid
39      !!       
40      !! ** Action  : - Initialisation of some variables
41      !!              - Definition of some constants linked with the grid
42      !!              - Definition of the metric coef. for the sea/ice
43      !!              - Initialization of the ice masks (tmsk, umsk)
44      !!
45      !! References : Deleersnijder et al. Ocean Modelling 100, 7-10
46      !!---------------------------------------------------------------------
47      INTEGER :: ji, jj      ! dummy loop indices
48      REAL(wp) ::   zusden   ! local scalars
49#if defined key_lim2_vp
50      REAL(wp) ::   zusden2           ! local scalars
51      REAL(wp) ::   zh1p  , zh2p      !   -      -
52      REAL(wp) ::   zd2d1p, zd1d2p    !   -      -
53      REAL(wp), DIMENSION(jpi,jpj) ::   zd2d1 , zd1d2   ! 2D workspace
54#endif
55      !!---------------------------------------------------------------------
56      !
57      IF(lwp) THEN
58         WRITE(numout,*)
59         WRITE(numout,*) 'lim_msh_2 : LIM 2.0 sea-ice model, mesh initialization'
60         WRITE(numout,*) '~~~~~~~~~'
61      ENDIF
62     
63      !----------------------------------------------------------                         
64      !    Initialization of local and some global (common) variables
65      !------------------------------------------------------------------
66     
67      njeq   = INT( jpj / 2 )   !i bug mpp potentiel
68      njeqm1 = njeq - 1 
69
70      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor at T-point
71 
72      IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN   ! local domain include both hemisphere
73         l_jeq = .TRUE.
74         njeq  = 1
75         DO WHILE ( njeq <= jpj .AND. fcor(1,njeq) < 0.e0 )
76            njeq = njeq + 1
77         END DO
78         IF(lwp ) WRITE(numout,*) '          the equator is inside the domain at about njeq = ', njeq
79      ELSEIF( fcor(1,1) < 0.e0 ) THEN
80         l_jeq = .FALSE.
81         njeq = jpj
82         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the southern hemisphere: njeq = ', njeq
83      ELSE
84         l_jeq = .FALSE.
85         njeq = 2
86         IF(lwp ) WRITE(numout,*) '          the model domain is entirely in the northern hemisphere: njeq = ', njeq
87      ENDIF
88
89      njeqm1 = njeq - 1
90
91
92      !   For each grid, definition of geometric tables
93      !------------------------------------------------------------------
94     
95      !-------------------
96      ! Conventions :    !
97      !-------------------
98      !  indices 1 \ 2 <-> localisation in the 2 direction x \ y
99      !  3rd indice <-> localisation on the mesh :
100      !  0 = Centre ;  1 = corner W x(i-1/2) ; 2 = corner S y(j-1/2) ;
101      !  3 = corner SW x(i-1/2),y(j-1/2)
102      !-------------------
103!!ibug ???
104      wght(:,:,:,:) = 0.e0
105      tmu(:,:)      = 0.e0
106#if defined key_lim2_vp 
107      akappa(:,:,:,:)     = 0.e0
108      alambd(:,:,:,:,:,:) = 0.e0
109#else
110      tmv(:,:) = 0.e0
111      tmf(:,:) = 0.e0
112#endif
113!!i
114     
115#if defined key_lim2_vp
116      ! metric coefficients for sea ice dynamic (VP rheology)
117      !----------------------------------------
118      !                                                       ! akappa
119      DO jj = 2, jpj
120         zd1d2(:,jj) = e1v(:,jj) - e1v(:,jj-1)
121      END DO
122      DO ji = 2, jpi
123         zd2d1(ji,:) = e2u(ji,:) - e2u(ji-1,:)
124      END DO
125      CALL lbc_lnk( zd1d2, 'T', -1. )   ;   CALL lbc_lnk( zd2d1, 'T', -1. )      ! lateral boundary condition
126      !
127      akappa(:,:,1,1) =        1.0 / ( 2.0 * e1t(:,:) )
128      akappa(:,:,1,2) = zd1d2(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) )
129      akappa(:,:,2,1) = zd2d1(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) )
130      akappa(:,:,2,2) =        1.0 / ( 2.0 * e2t(:,:) )
131     
132      !             
133      DO jj = 2, jpj                                         ! weights (wght) at I-points
134         DO ji = 2, jpi
135            zusden = 1. / (  ( e1t(ji,jj) + e1t(ji-1,jj  ) )   &
136               &           * ( e2t(ji,jj) + e2t(ji  ,jj-1) ) )
137            wght(ji,jj,1,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj  )
138            wght(ji,jj,1,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj-1)
139            wght(ji,jj,2,1) = zusden * e1t(ji-1,jj) * e2t(ji,jj  )
140            wght(ji,jj,2,2) = zusden * e1t(ji-1,jj) * e2t(ji,jj-1)
141         END DO
142      END DO
143      CALL lbc_lnk( wght(:,:,1,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V-point
144      CALL lbc_lnk( wght(:,:,1,2), 'I', 1. )      ! the value of wght at jpj is wrong
145      CALL lbc_lnk( wght(:,:,2,1), 'I', 1. )      ! but it is never used
146      CALL lbc_lnk( wght(:,:,2,2), 'I', 1. )
147#else
148      ! metric coefficients for sea ice dynamic (EVP rheology)
149      !----------------------------------------
150      DO jj = 1, jpjm1                                       ! weights (wght) at F-points
151         DO ji = 1, jpim1
152            zusden = 1. / (  ( e1t(ji+1,jj  ) + e1t(ji,jj) )   &
153               &           * ( e2t(ji  ,jj+1) + e2t(ji,jj) ) ) 
154            wght(ji,jj,1,1) = zusden * e1t(ji+1,jj) * e2t(ji,jj+1)
155            wght(ji,jj,1,2) = zusden * e1t(ji+1,jj) * e2t(ji,jj  )
156            wght(ji,jj,2,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj+1)
157            wght(ji,jj,2,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj  )
158         END DO
159      END DO
160      CALL lbc_lnk( wght(:,:,1,1), 'F', 1. )   ;   CALL lbc_lnk( wght(:,:,1,2), 'F', 1. )       ! lateral boundary cond.   
161      CALL lbc_lnk( wght(:,:,2,1), 'F', 1. )   ;   CALL lbc_lnk( wght(:,:,2,2), 'F', 1. )
162#endif
163   
164      ! Coefficients for divergence of the stress tensor
165      !-------------------------------------------------
166#if defined key_lim2_vp
167      DO jj = 2, jpj                                     ! VP rheology
168         DO ji = 2, jpi   ! NO vector opt.
169            zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2) + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   &
170               &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1) + e1t(ji-1,jj-1) * wght(ji,jj,1,1)
171            zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2) + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   &
172               &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1) + e2t(ji-1,jj-1) * wght(ji,jj,1,1)
173            !
174            zusden  = 1.e0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 )
175            zusden2 = zusden * 2.0 
176            zd1d2p  = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   )
177            zd2d1p  = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   )
178            !
179            alambd(ji,jj,2,2,2,1) = zusden2 * e2t(ji  ,jj-1)
180            alambd(ji,jj,2,2,2,2) = zusden2 * e2t(ji  ,jj  )
181            alambd(ji,jj,2,2,1,1) = zusden2 * e2t(ji-1,jj-1)
182            alambd(ji,jj,2,2,1,2) = zusden2 * e2t(ji-1,jj  )
183            !
184            alambd(ji,jj,1,1,2,1) = zusden2 * e1t(ji  ,jj-1)
185            alambd(ji,jj,1,1,2,2) = zusden2 * e1t(ji  ,jj  )
186            alambd(ji,jj,1,1,1,1) = zusden2 * e1t(ji-1,jj-1)
187            alambd(ji,jj,1,1,1,2) = zusden2 * e1t(ji-1,jj  )
188            !
189            alambd(ji,jj,1,2,2,1) = zd1d2p
190            alambd(ji,jj,1,2,2,2) = zd1d2p
191            alambd(ji,jj,1,2,1,1) = zd1d2p
192            alambd(ji,jj,1,2,1,2) = zd1d2p
193            !
194            alambd(ji,jj,2,1,2,1) = zd2d1p
195            alambd(ji,jj,2,1,2,2) = zd2d1p
196            alambd(ji,jj,2,1,1,1) = zd2d1p
197            alambd(ji,jj,2,1,1,2) = zd2d1p
198         END DO
199      END DO
200      !
201      ! lateral boundary conditions
202      ! CAUTION: even with the lbc_lnk at ice U-V point, the value of wght at jpj is wrong but it is never used
203      CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. ) 
204      CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. ) 
205      !
206      CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. )
207      CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. )
208      !
209      CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. )
210      CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. )
211      !
212      CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. )
213      CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. )
214#endif
215           
216
217      ! Initialization of ice masks
218      !----------------------------
219      !     
220      tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask
221      !
222#if defined key_lim2_vp
223      ! VP rheology : ice velocity point is I-point
224      tmu(:,1) = 0.e0              !
225      tmu(1,:) = 0.e0
226      DO jj = 2, jpj
227         DO ji = 2, jpim1   ! NO vector opt.
228            tmu(ji,jj) =  tms(ji,jj) * tms(ji-1,jj) * tms(ji,jj-1) * tms(ji-1,jj-1)           
229         END DO
230      END DO
231      CALL lbc_lnk( tmu(:,:), 'I', 1. )      ! lateral boundary conditions   
232#else
233      ! EVP rheology : ice velocity point are U- & V-points ; ice vorticity point is F-point
234      tmu(:,:) = umask(:,:,1)
235      tmv(:,:) = vmask(:,:,1)
236      tmf(:,:) = 0.e0                        ! used of fmask except its special value along the coast (rn_shlat)
237      WHERE( fmask(:,:,1) == 1.e0 )   tmf(:,:) = 1.e0
238#endif
239      !
240      area(:,:) = e1t(:,:) * e2t(:,:)        ! unmasked and masked area of T-grid cell
241      !
242   END SUBROUTINE lim_msh_2
243
244#else
245   !!----------------------------------------------------------------------
246   !!   Default option            Dummy Module         NO LIM sea-ice model
247   !!----------------------------------------------------------------------
248CONTAINS
249   SUBROUTINE lim_msh_2           ! Dummy routine
250   END SUBROUTINE lim_msh_2
251#endif
252
253   !!======================================================================
254END MODULE limmsh_2
Note: See TracBrowser for help on using the repository browser.