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/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

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