source: branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/C1D/domc1d.F90 @ 4245

Last change on this file since 4245 was 4245, checked in by cetlod, 8 years ago

dev_locean_cmcc_ingv_ukmo_merc : merge in the MERC_UKMO dev branch with trunk rev 4119

File size: 8.9 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 iom                           ! I/O library (iom_get)
16   USE in_out_manager                ! I/O manager (ctmp1)
17   USE dom_oce , ONLY : nimpp, njmpp ! Shared/distributed memory setting (mpp_init routine)
18   USE wrk_nemo                      ! Memory allocation
19   USE timing                        ! Timing
20
21   IMPLICIT NONE
22   PRIVATE
23
24   PUBLIC   dom_c1d                  ! Routine called in domcfg.F90
25
26   !!----------------------------------------------------------------------
27   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
28   !! $Id: domc1d.F90 3851 2013-04-30 10:30:51Z hadcv $
29   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
30   !!----------------------------------------------------------------------
31CONTAINS
32   
33   SUBROUTINE dom_c1d( plat, plon )
34      !!----------------------------------------------------------------------
35      !!                   ***  ROUTINE dom_c1d  ***
36      !!
37      !! ** Purpose : Recalculate jpizoom/jpjzoom indices from lat/lon point
38      !!
39      !! ** Method  : Calculate global gphit and glamt as for dom_hgr.
40      !!              After, find closest grid point to lat/lon point as for
41      !!              dom_ngb on T grid. From this infer jpizoom and jpjzoom.
42      !!
43      !! ** Action  : Recalculate jpizoom, jpjzoom (indices of C1D zoom)
44      !!----------------------------------------------------------------------
45      INTEGER  ::  ji, jj                          ! Dummy loop indices
46      INTEGER  ::  inum                            ! Coordinate file handle (case 0)
47      INTEGER  ::  ijeq                            ! Index of equator T point (case 4)
48
49      INTEGER , DIMENSION(2) ::   iloc             ! Minloc returned indices
50
51      REAL(wp), INTENT(in) ::  plat                ! Column latitude
52      REAL(wp), INTENT(in) ::  plon                ! Column longitude
53
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
60      REAL(wp) , POINTER, DIMENSION(:,:) ::  gphidta, glamdta, zdist ! Global lat/lon
61      !!----------------------------------------------------------------------
62
63      IF( nn_timing == 1 )   CALL timing_start('dom_c1d')
64
65      CALL wrk_alloc( jpidta, jpjdta, gphidta, glamdta, zdist )
66
67
68      ! ============================= !
69      !  Code from dom_hgr:           !
70      !  Calculate global horizontal  !
71      !  mesh, only glamt and gphit   !
72      ! ============================= !
73
74      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
75
76      CASE ( 0 )                 !  curvilinear coordinate on the sphere read in coordinate.nc file
77
78         CALL iom_open( 'coordinates', inum )
79         CALL iom_get( inum, jpdom_unknown, 'glamt', glamdta ) ! mig, mjg undefined at this point
80         CALL iom_get( inum, jpdom_unknown, 'gphit', gphidta ) ! so use jpdom_unknown not jpdom_data
81         CALL iom_close ( inum )
82
83         PRINT *,'Check dom_c1d coordinates file data read in:' !!!
84         PRINT *,'Bottom-left most glamdta is ', glamdta(1,1)    !!! Need to check
85         PRINT *,'Bottom-left most gphidta is ', gphidta(1,1)    !!! field read
86         PRINT *,'We are using nimpp,njmpp = ' , nimpp,njmpp     !!!
87
88      CASE ( 1 )                 ! geographical mesh on the sphere with regular grid-spacing
89
90         DO jj = 1, jpjdta
91            DO ji = 1, jpidta
92               zti = FLOAT( ji - 1 + nimpp - 1 )
93               ztj = FLOAT( jj - 1 + njmpp - 1 )
94
95               glamdta(ji,jj) = ppglam0 + ppe1_deg * zti
96               gphidta(ji,jj) = ppgphi0 + ppe2_deg * ztj
97            END DO
98         END DO
99
100      CASE ( 2:3 )               ! f- or beta-plane with regular grid-spacing
101         
102         glam0 = 0.e0
103         gphi0 = - ppe2_m * 1.e-3
104
105         DO jj = 1, jpjdta
106            DO ji = 1, jpidta
107               glamdta(ji,jj) = glam0 + ppe1_m * 1.e-3 * FLOAT( ji - 1 + nimpp - 1 )
108               gphidta(ji,jj) = gphi0 + ppe2_m * 1.e-3 * FLOAT( jj - 1 + njmpp - 1 )
109            END DO
110         END DO
111
112      CASE ( 4 )                 ! geographical mesh on the sphere, isotropic MERCATOR type
113
114         IF( ppgphi0 == -90 )   CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
115
116         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
117         ijeq = ABS( 180. / rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
118         IF( ppgphi0 > 0 )   ijeq = -ijeq
119
120         DO jj = 1, jpjdta
121            DO ji = 1, jpidta
122               zti = FLOAT( ji - 1    + nimpp - 1 )
123               ztj = FLOAT( jj - ijeq + njmpp - 1 )
124
125               glamdta(ji,jj) = ppglam0 + ppe1_deg * zti
126               gphidta(ji,jj) = 1. / rad * ASIN ( TANH( ppe1_deg * rad * ztj ) )
127            END DO
128         END DO
129
130      CASE ( 5 )                 ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration)
131   
132         zlam1 = -85
133         zphi1 = 29
134         ze1 = 106000. / FLOAT(jp_cfg)
135 
136         zsin_alpha = - SQRT( 2. ) / 2.
137         zcos_alpha =   SQRT( 2. ) / 2.
138         ze1deg = ze1 / (ra * rad)
139
140         glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjdta-2 ) ! Force global
141         gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjdta-2 )
142
143         DO jj = 1, jpjdta
144            DO ji = 1, jpidta
145               zim05 = FLOAT( ji + nimpp - 1 ) - 1.5
146               zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5
147
148               glamdta(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
149               gphidta(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
150            END DO
151         END DO
152
153      CASE DEFAULT
154
155         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
156         CALL ctl_stop( ctmp1 )
157
158      END SELECT
159
160
161      ! ============================== !
162      !  Code from dom_ngb:            !
163      !  Calculate the nearest grid    !
164      !  point to the given lat/lon &  !
165      !  update jpizoom and jpjzoom    !
166      ! ============================== !
167
168      zlon         = MOD( plon         + 720., 360. )                                      ! plon    between    0 and 360
169      glamdta(:,:) = MOD( glamdta(:,:) + 720., 360. )                                      ! glamdta between    0 and 360
170      IF( zlon > 270. )   zlon = zlon - 360.                                               ! zlon    between  -90 and 270
171      IF( zlon <  90. )   WHERE( glamdta(:,:) > 180. ) glamdta(:,:) = glamdta(:,:) - 360.  ! glamdta between -180 and 180
172
173      glamdta(:,:) = glamdta(:,:) - zlon
174      gphidta(:,:) = gphidta(:,:) - plat
175      zdist(:,:)   = glamdta(:,:) * glamdta(:,:) + gphidta(:,:) * gphidta(:,:)
176     
177      iloc(:) = MINLOC( zdist(:,:) ) ! No mask; zoom indices freely defined
178      jpizoom = iloc(1) + nimpp - 2  ! Minloc index - 1; want the bottom-left
179      jpjzoom = iloc(2) + njmpp - 2  ! corner index of the zoom domain.
180
181      CALL wrk_dealloc( jpidta, jpjdta, gphidta, glamdta, zdist )
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      IF( nn_timing == 1 )   CALL timing_stop('dom_c1d')
193
194   END SUBROUTINE dom_c1d
195
196#else
197   !!----------------------------------------------------------------------
198   !!   Default option                                  NO 1D Configuration
199   !!----------------------------------------------------------------------
200CONTAINS 
201   SUBROUTINE dom_c1d( plat, plon )     ! Empty routine
202      REAL :: plat, plon
203      WRITE(*,*) 'dom_c1d: You should not have seen this print! error?',plat,plon
204   END SUBROUTINE dom_c1d
205#endif
206
207   !!======================================================================
208END MODULE domc1d
Note: See TracBrowser for help on using the repository browser.