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.
domutl.F90 in NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DOM – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DOM/domutl.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 3 years ago

Mixed precision version, tested up to 30 years on ORCA2.
  • Property svn:keywords set to Id
File size: 7.5 KB
Line 
1MODULE domutl
2   !!======================================================================
3   !!                       ***  MODULE  domutl  ***
4   !! Grid utilities:
5   !!======================================================================
6   !! History : 4.2  !  2020-04  (S. Masson)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   dom_ngb       : find the closest grid point from a given lon/lat position
11   !!   dom_uniq      : identify unique point of a grid (TUVF)
12   !!----------------------------------------------------------------------
13   !
14   USE dom_oce        ! ocean space and time domain
15   !
16   USE in_out_manager ! I/O manager
17   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
18   USE lib_mpp        ! for mppsum
19
20   IMPLICIT NONE
21   PRIVATE
22
23   INTERFACE is_tile
24      MODULE PROCEDURE is_tile_2d_sp, is_tile_3d_sp, is_tile_4d_sp, is_tile_2d_dp, is_tile_3d_dp, is_tile_4d_dp
25   END INTERFACE is_tile
26
27   PUBLIC dom_ngb    ! routine called in iom.F90 module
28   PUBLIC dom_uniq   ! Called by dommsk and domwri
29   PUBLIC is_tile
30
31   !!----------------------------------------------------------------------
32   !! NEMO/OCE 4.2 , NEMO Consortium (2020)
33   !! $Id$
34   !! Software governed by the CeCILL license (see ./LICENSE)
35   !!----------------------------------------------------------------------
36#  include "single_precision_substitute.h90"
37CONTAINS
38
39   SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, kkk )
40      !!----------------------------------------------------------------------
41      !!                    ***  ROUTINE dom_ngb  ***
42      !!
43      !! ** Purpose :   find the closest grid point from a given lon/lat position
44      !!
45      !! ** Method  :   look for minimum distance in cylindrical projection
46      !!                -> not good if located at too high latitude...
47      !!----------------------------------------------------------------------
48      REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point
49      INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point
50      INTEGER         , INTENT(in   ), OPTIONAL :: kkk  ! k-index of the mask level used
51      CHARACTER(len=1), INTENT(in   ) ::   cdgrid       ! grid name 'T', 'U', 'V', 'W'
52      !
53      INTEGER :: ik         ! working level
54      INTEGER , DIMENSION(2) ::   iloc
55      REAL(wp)                :: zlon
56      REAL(dp)                :: zmini
57      REAL(wp), DIMENSION(jpi,jpj) ::   zglam, zgphi, zdist
58      LOGICAL , DIMENSION(jpi,jpj) ::   llmsk
59      !!--------------------------------------------------------------------
60      !
61      ik = 1
62      IF ( PRESENT(kkk) ) ik=kkk
63      !
64      SELECT CASE( cdgrid )
65      CASE( 'U' ) ;   zglam(:,:) = glamu(:,:)   ;   zgphi(:,:) = gphiu(:,:)   ;   llmsk(:,:) = tmask_h(:,:) * umask(:,:,ik) == 1._wp
66      CASE( 'V' ) ;   zglam(:,:) = glamv(:,:)   ;   zgphi(:,:) = gphiv(:,:)   ;   llmsk(:,:) = tmask_h(:,:) * vmask(:,:,ik) == 1._wp
67      CASE( 'F' ) ;   zglam(:,:) = glamf(:,:)   ;   zgphi(:,:) = gphif(:,:)   ;   llmsk(:,:) = tmask_h(:,:) * fmask(:,:,ik) == 1._wp
68      CASE DEFAULT;   zglam(:,:) = glamt(:,:)   ;   zgphi(:,:) = gphit(:,:)   ;   llmsk(:,:) = tmask_h(:,:) * tmask(:,:,ik) == 1._wp
69      END SELECT
70      !
71      zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360
72      zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360
73      IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270
74      IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180
75      zglam(:,:) = zglam(:,:) - zlon
76      !
77      zgphi(:,:) = zgphi(:,:) - plat
78      zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:)
79      !
80      CALL mpp_minloc( 'domngb', CASTDP(zdist(:,:)), llmsk, zmini, iloc, ldhalo = .TRUE. )
81      kii = iloc(1)
82      kjj = iloc(2)
83      !
84   END SUBROUTINE dom_ngb
85
86
87   SUBROUTINE dom_uniq( puniq, cdgrd )
88      !!----------------------------------------------------------------------
89      !!                  ***  ROUTINE dom_uniq  ***
90      !!
91      !! ** Purpose :   identify unique point of a grid (TUVF)
92      !!
93      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element
94      !!                2) check which elements have been changed
95      !!----------------------------------------------------------------------
96      CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   !
97      REAL(dp), DIMENSION(:,:), INTENT(inout) ::   puniq   !
98      !
99      REAL(wp)                       ::  zshift   ! shift value link to the process number
100      INTEGER                        ::  ji       ! dummy loop indices
101      LOGICAL , DIMENSION(jpi,jpj,1) ::   lluniq  ! store whether each point is unique or not
102      REAL(wp), DIMENSION(jpi,jpj  ) ::   ztstref
103      !!----------------------------------------------------------------------
104      !
105      ! build an array with different values for each element
106      ! in mpp: make sure that these values are different even between process
107      ! -> apply a shift value according to the process number
108      zshift = jpimax * jpjmax * ( narea - 1 )
109      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) )
110      !
111      puniq(:,:) = ztstref(:,:)                    ! default definition
112      CALL lbc_lnk( 'domwri', puniq, cdgrd, 1._wp )   ! apply boundary conditions
113      lluniq(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have not been changed
114      !
115      puniq(:,:) = REAL( COUNT( lluniq(:,:,:), dim = 3 ), wp )
116      !
117   END SUBROUTINE dom_uniq
118
119
120   INTEGER FUNCTION is_tile_2d_sp( pt )
121      REAL(sp), DIMENSION(:,:), INTENT(in) ::   pt
122
123      IF( l_istiled .AND. (SIZE(pt, 1) < jpi .OR. SIZE(pt, 2) < jpj) ) THEN
124         is_tile_2d_sp = 1
125      ELSE
126         is_tile_2d_sp = 0
127      ENDIF
128   END FUNCTION is_tile_2d_sp
129
130
131   INTEGER FUNCTION is_tile_2d_dp( pt )
132      REAL(dp), DIMENSION(:,:), INTENT(in) ::   pt
133
134      IF( l_istiled .AND. (SIZE(pt, 1) < jpi .OR. SIZE(pt, 2) < jpj) ) THEN
135         is_tile_2d_dp = 1
136      ELSE
137         is_tile_2d_dp = 0
138      ENDIF
139   END FUNCTION is_tile_2d_dp
140
141
142   INTEGER FUNCTION is_tile_3d_sp( pt )
143      REAL(sp), DIMENSION(:,:,:), INTENT(in) ::   pt
144
145      IF( l_istiled .AND. (SIZE(pt, 1) < jpi .OR. SIZE(pt, 2) < jpj) ) THEN
146         is_tile_3d_sp = 1
147      ELSE
148         is_tile_3d_sp = 0
149      ENDIF
150   END FUNCTION is_tile_3d_sp
151
152
153   INTEGER FUNCTION is_tile_3d_dp( pt )
154      REAL(dp), DIMENSION(:,:,:), INTENT(in) ::   pt
155
156      IF( l_istiled .AND. (SIZE(pt, 1) < jpi .OR. SIZE(pt, 2) < jpj) ) THEN
157         is_tile_3d_dp = 1
158      ELSE
159         is_tile_3d_dp = 0
160      ENDIF
161   END FUNCTION is_tile_3d_dp
162
163
164   INTEGER FUNCTION is_tile_4d_sp( pt )
165      REAL(sp), DIMENSION(:,:,:,:), INTENT(in) ::   pt
166
167      IF( l_istiled .AND. (SIZE(pt, 1) < jpi .OR. SIZE(pt, 2) < jpj) ) THEN
168         is_tile_4d_sp = 1
169      ELSE
170         is_tile_4d_sp = 0
171      ENDIF
172   END FUNCTION is_tile_4d_sp
173
174
175   INTEGER FUNCTION is_tile_4d_dp( pt )
176      REAL(dp), DIMENSION(:,:,:,:), INTENT(in) ::   pt
177
178      IF( l_istiled .AND. (SIZE(pt, 1) < jpi .OR. SIZE(pt, 2) < jpj) ) THEN
179         is_tile_4d_dp = 1
180      ELSE
181         is_tile_4d_dp = 0
182      ENDIF
183   END FUNCTION is_tile_4d_dp
184   !!======================================================================
185END MODULE domutl
Note: See TracBrowser for help on using the repository browser.