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/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90 @ 7512

Last change on this file since 7512 was 7508, checked in by mocavero, 8 years ago

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