source: NEMO/trunk/src/OCE/C1D/domc1d.F90 @ 9598

Last change on this file since 9598 was 9598, checked in by nicolasmartin, 2 years ago

Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure
Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing
Standardisation of routines header with version 4.0 and year 2018
Fix for some broken symlinks

  • Property svn:keywords set to Id
File size: 9.3 KB
Line 
1MODULE domc1d
2   !!======================================================================
3   !!                     ***  MODULE  domc1d  ***
4   !! Ocean Domain : 1D column position from lat/lon namelist specification
5   !!======================================================================
6   !! History :  3.5  !  2013-04  (D. Calvert)  Original code
7   !!----------------------------------------------------------------------
8#if defined key_c1d
9   !!----------------------------------------------------------------------
10   !!   'key_c1d'   :                                      1D Configuration
11   !!----------------------------------------------------------------------
12   !!   dom_c1d     : Determine jpizoom/jpjzoom from a given lat/lon
13   !!----------------------------------------------------------------------
14   USE phycst         ! Physical constants (and par_oce)
15   USE dom_oce , ONLY : nimpp, njmpp ! Shared/distributed memory setting
16   !
17   USE iom            ! I/O library (iom_get)
18   USE in_out_manager ! I/O manager (ctmp1)
19
20   IMPLICIT NONE
21   PRIVATE
22
23   PUBLIC   dom_c1d   ! called in domcfg.F90
24
25   INTEGER ::   jpizoom = 1      !: left bottom (i,j) indices of the zoom
26   INTEGER ::   jpjzoom = 1      !: in data domain indices
27
28   !!----------------------------------------------------------------------
29   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
30   !! $Id$
31   !! Software governed by the CeCILL licence     (./LICENSE)
32   !!----------------------------------------------------------------------
33CONTAINS
34   
35   SUBROUTINE dom_c1d( plat, plon )
36      !!----------------------------------------------------------------------
37      !!                   ***  ROUTINE dom_c1d  ***
38      !!
39      !! ** Purpose : Recalculate jpizoom/jpjzoom indices from lat/lon point
40      !!
41      !! ** Method  : Calculate global gphit and glamt as for dom_hgr.
42      !!              After, find closest grid point to lat/lon point as for
43      !!              dom_ngb on T grid. From this infer jpizoom and jpjzoom.
44      !!
45      !! ** Action  : Recalculate jpizoom, jpjzoom (indices of C1D zoom)
46      !!----------------------------------------------------------------------
47      REAL(wp), INTENT(in) ::  plat, plon    ! Column latitude &  longitude
48      !
49      INTEGER  ::  ji, jj   ! Dummy loop indices
50      INTEGER  ::  inum     ! Coordinate file handle (case 0)
51      INTEGER  ::  ijeq     ! Index of equator T point (case 4)
52      INTEGER  ::  ios      ! Local integer output status for namelist read
53      INTEGER , DIMENSION(2) ::   iloc   ! Minloc returned indices
54      REAL(wp) ::  zlon                            ! Wraparound longitude
55      REAL(wp) ::  zti, ztj, zarg                  ! Local scalars
56      REAL(wp) ::  glam0, gphi0                    ! Variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
57      REAL(wp) ::  zlam1, zcos_alpha, ze1, ze1deg  ! Case 5 local scalars
58      REAL(wp) ::  zphi1, zsin_alpha, zim05, zjm05 !         
59      REAL(wp) , DIMENSION(jpidta,jpjdta) ::  gphidta, glamdta, zdist ! Global lat/lon
60      !!
61      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, ln_meshmask, rn_hmin,   &
62         &             rn_atfp     , rn_rdt , ln_crs,  jphgr_msh, &
63         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
64         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
65         &             ppa2, ppkth2, ppacr2
66      !!----------------------------------------------------------------------
67
68      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
69      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 901 )
70901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
71      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
72      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 902 )
73902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
74
75
76      ! ============================= !
77      !  Code from dom_hgr:           !
78      !  Calculate global horizontal  !
79      !  mesh, only glamt and gphit   !
80      ! ============================= !
81      !
82      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
83      !
84      CASE ( 0 )                 !  curvilinear coordinate on the sphere read in coordinate.nc file
85         !
86         CALL iom_open( 'coordinates', inum )
87         CALL iom_get( inum, jpdom_unknown, 'glamt', glamdta ) ! mig, mjg undefined at this point
88         CALL iom_get( inum, jpdom_unknown, 'gphit', gphidta ) ! so use jpdom_unknown not jpdom_data
89         CALL iom_close ( inum )
90         !
91      CASE ( 1 )                 ! geographical mesh on the sphere with regular grid-spacing
92         !
93         DO jj = 1, jpjdta
94            DO ji = 1, jpidta
95               zti = FLOAT( ji - 1 + nimpp - 1 )
96               ztj = FLOAT( jj - 1 + njmpp - 1 )
97               !
98               glamdta(ji,jj) = ppglam0 + ppe1_deg * zti
99               gphidta(ji,jj) = ppgphi0 + ppe2_deg * ztj
100            END DO
101         END DO
102         !
103      CASE ( 2:3 )               ! f- or beta-plane with regular grid-spacing
104         !
105         glam0 = 0.e0
106         gphi0 = - ppe2_m * 1.e-3
107         !
108         DO jj = 1, jpjdta
109            DO ji = 1, jpidta
110               glamdta(ji,jj) = glam0 + ppe1_m * 1.e-3 * FLOAT( ji - 1 + nimpp - 1 )
111               gphidta(ji,jj) = gphi0 + ppe2_m * 1.e-3 * FLOAT( jj - 1 + njmpp - 1 )
112            END DO
113         END DO
114         !
115      CASE ( 4 )                 ! geographical mesh on the sphere, isotropic MERCATOR type
116         !
117         IF( ppgphi0 == -90 )   CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
118         !
119         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
120         ijeq = ABS( 180. / rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
121         IF( ppgphi0 > 0 )   ijeq = -ijeq
122         !
123         DO jj = 1, jpjdta
124            DO ji = 1, jpidta
125               zti = FLOAT( ji - 1    + nimpp - 1 )
126               ztj = FLOAT( jj - ijeq + njmpp - 1 )
127               !
128               glamdta(ji,jj) = ppglam0 + ppe1_deg * zti
129               gphidta(ji,jj) = 1. / rad * ASIN ( TANH( ppe1_deg * rad * ztj ) )
130            END DO
131         END DO
132         !
133      CASE ( 5 )                 ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
134         !
135         zlam1 = -85
136         zphi1 = 29
137         ze1 = 106000. / REAL( nn_cfg , wp )
138         !
139         zsin_alpha = - SQRT( 2. ) / 2.
140         zcos_alpha =   SQRT( 2. ) / 2.
141         ze1deg = ze1 / (ra * rad)
142         !
143         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjdta-2 ) ! Force global
144         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjdta-2 )
145         !
146         DO jj = 1, jpjdta
147            DO ji = 1, jpidta
148               zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
149               zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
150
151               glamdta(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
152               gphidta(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
153            END DO
154         END DO
155         !
156      CASE DEFAULT
157         !
158         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
159         CALL ctl_stop( ctmp1 )
160         !
161      END SELECT
162
163      ! ============================== !
164      !  Code from dom_ngb:            !
165      !  Calculate the nearest grid    !
166      !  point to the given lat/lon &  !
167      !  update jpizoom and jpjzoom    !
168      ! ============================== !
169
170      zlon         = MOD( plon         + 720., 360. )                                      ! plon    between    0 and 360
171      glamdta(:,:) = MOD( glamdta(:,:) + 720., 360. )                                      ! glamdta between    0 and 360
172      IF( zlon > 270. )   zlon = zlon - 360.                                               ! zlon    between  -90 and 270
173      IF( zlon <  90. )   WHERE( glamdta(:,:) > 180. ) glamdta(:,:) = glamdta(:,:) - 360.  ! glamdta between -180 and 180
174
175      glamdta(:,:) = glamdta(:,:) - zlon
176      gphidta(:,:) = gphidta(:,:) - plat
177      zdist(:,:)   = glamdta(:,:) * glamdta(:,:) + gphidta(:,:) * gphidta(:,:)
178     
179      iloc(:) = MINLOC( zdist(:,:) ) ! No mask; zoom indices freely defined
180      jpizoom = iloc(1) + nimpp - 2  ! Minloc index - 1; want the bottom-left
181      jpjzoom = iloc(2) + njmpp - 2  ! corner index of the zoom domain.
182
183      IF(lwp) THEN
184         WRITE(numout,*)
185         WRITE(numout,*) 'dom_c1d : compute jpizoom & jpjzoom from global mesh and given coordinates'
186         WRITE(numout,*) '~~~~~~~'
187         WRITE(numout,*) '      column i zoom index             jpizoom = ', jpizoom
188         WRITE(numout,*) '      column j zoom index             jpjzoom = ', jpjzoom
189         WRITE(numout,*)
190      ENDIF
191      !
192   END SUBROUTINE dom_c1d
193
194#else
195   !!----------------------------------------------------------------------
196   !!   Default option                                  NO 1D Configuration
197   !!----------------------------------------------------------------------
198CONTAINS 
199   SUBROUTINE dom_c1d( plat, plon )     ! Empty routine
200      REAL :: plat, plon
201      WRITE(*,*) 'dom_c1d: You should not have seen this print! error?',plat,plon
202   END SUBROUTINE dom_c1d
203#endif
204
205   !!======================================================================
206END MODULE domc1d
Note: See TracBrowser for help on using the repository browser.