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.
geo2ocean.F90 in branches/2013/dev_r3987_UKMO4_OBS/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2013/dev_r3987_UKMO4_OBS/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 @ 4113

Last change on this file since 4113 was 4113, checked in by andrewryan, 11 years ago

Merged with trunk at revision 4103 to keep branch up to date and tested with latest developments
  • Property svn:keywords set to Id
File size: 26.2 KB
Line 
1MODULE geo2ocean
2   !!======================================================================
3   !!                     ***  MODULE  geo2ocean  ***
4   !! Ocean mesh    :  ???
5   !!======================================================================
6   !! History :  OPA  !  07-1996  (O. Marti)  Original code
7   !!   NEMO     1.0  !  02-2008  (G. Madec)  F90: Free form
8   !!            3.0  ! 
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   repcmo      :
13   !!   angle       :
14   !!   geo2oce     :
15   !!   repere      :   old routine suppress it ???
16   !!----------------------------------------------------------------------
17   USE dom_oce         ! mesh and scale factors
18   USE phycst          ! physical constants
19   USE in_out_manager  ! I/O manager
20   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
21   USE lib_mpp         ! MPP library
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   rot_rep, repcmo, repere, geo2oce, oce2geo   ! only rot_rep should be used
27                                             ! repcmo and repere are keep only for compatibility.
28                                             ! they are only a useless overlay of rot_rep
29
30   PUBLIC   obs_rot
31
32   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   &
33      gsint, gcost,   &  ! cos/sin between model grid lines and NP direction at T point
34      gsinu, gcosu,   &  ! cos/sin between model grid lines and NP direction at U point
35      gsinv, gcosv,   &  ! cos/sin between model grid lines and NP direction at V point
36      gsinf, gcosf       ! cos/sin between model grid lines and NP direction at F point
37
38   LOGICAL ,              SAVE, DIMENSION(4)     ::   linit = .FALSE.
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gsinlon, gcoslon, gsinlat, gcoslat
40
41   LOGICAL ::   lmust_init = .TRUE.        !: used to initialize the cos/sin variables (se above)
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1,   &
53                       px2 , py2 )
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE repcmo  ***
56      !!
57      !! ** Purpose :   Change vector componantes from a geographic grid to a
58      !!      stretched coordinates grid.
59      !!
60      !! ** Method  :   Initialization of arrays at the first call.
61      !!
62      !! ** Action  : - px2 : first  componante (defined at u point)
63      !!              - py2 : second componante (defined at v point)
64      !!----------------------------------------------------------------------
65      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pxu1, pyu1   ! geographic vector componantes at u-point
66      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pxv1, pyv1   ! geographic vector componantes at v-point
67      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   px2          ! i-componante (defined at u-point)
68      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   py2          ! j-componante (defined at v-point)
69      !!----------------------------------------------------------------------
70     
71      ! Change from geographic to stretched coordinate
72      ! ----------------------------------------------
73      CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )
74      CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )
75     
76   END SUBROUTINE repcmo
77
78
79   SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot )
80      !!----------------------------------------------------------------------
81      !!                  ***  ROUTINE rot_rep  ***
82      !!
83      !! ** Purpose :   Rotate the Repere: Change vector componantes between
84      !!                geographic grid <--> stretched coordinates grid.
85      !!
86      !! History :
87      !!   9.2  !  07-04  (S. Masson) 
88      !!                  (O. Marti ) Original code (repere and repcmo)
89      !!----------------------------------------------------------------------
90      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) ::   pxin, pyin   ! vector componantes
91      CHARACTER(len=1),             INTENT( IN ) ::   cd_type      ! define the nature of pt2d array grid-points
92      CHARACTER(len=5),             INTENT( IN ) ::   cdtodo       ! specify the work to do:
93      !!                                                           ! 'en->i' east-north componantes to model i componante
94      !!                                                           ! 'en->j' east-north componantes to model j componante
95      !!                                                           ! 'ij->e' model i-j componantes to east componante
96      !!                                                           ! 'ij->n' model i-j componantes to east componante
97      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) ::   prot     
98      !!----------------------------------------------------------------------
99
100      ! Initialization of gsin* and gcos* at first call
101      ! -----------------------------------------------
102
103      IF( lmust_init ) THEN
104         IF(lwp) WRITE(numout,*)
105         IF(lwp) WRITE(numout,*) ' rot_rep : geographic <--> stretched'
106         IF(lwp) WRITE(numout,*) ' ~~~~~    coordinate transformation'
107         !
108         CALL angle       ! initialization of the transformation
109         lmust_init = .FALSE.
110      ENDIF
111     
112      SELECT CASE (cdtodo)
113      CASE ('en->i')      ! 'en->i' est-north componantes to model i componante
114         SELECT CASE (cd_type)
115         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:)
116         CASE ('U')   ;   prot(:,:) = pxin(:,:) * gcosu(:,:) + pyin(:,:) * gsinu(:,:)
117         CASE ('V')   ;   prot(:,:) = pxin(:,:) * gcosv(:,:) + pyin(:,:) * gsinv(:,:)
118         CASE ('F')   ;   prot(:,:) = pxin(:,:) * gcosf(:,:) + pyin(:,:) * gsinf(:,:)
119         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
120         END SELECT
121      CASE ('en->j')      ! 'en->j' est-north componantes to model j componante
122         SELECT CASE (cd_type)
123         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:)
124         CASE ('U')   ;   prot(:,:) = pyin(:,:) * gcosu(:,:) - pxin(:,:) * gsinu(:,:)
125         CASE ('V')   ;   prot(:,:) = pyin(:,:) * gcosv(:,:) - pxin(:,:) * gsinv(:,:)   
126         CASE ('F')   ;   prot(:,:) = pyin(:,:) * gcosf(:,:) - pxin(:,:) * gsinf(:,:)   
127         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
128         END SELECT
129      CASE ('ij->e')      ! 'ij->e' model i-j componantes to est componante
130         SELECT CASE (cd_type)
131         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:)
132         CASE ('U')   ;   prot(:,:) = pxin(:,:) * gcosu(:,:) - pyin(:,:) * gsinu(:,:)
133         CASE ('V')   ;   prot(:,:) = pxin(:,:) * gcosv(:,:) - pyin(:,:) * gsinv(:,:)
134         CASE ('F')   ;   prot(:,:) = pxin(:,:) * gcosf(:,:) - pyin(:,:) * gsinf(:,:)
135         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
136         END SELECT
137      CASE ('ij->n')      ! 'ij->n' model i-j componantes to est componante
138         SELECT CASE (cd_type)
139         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:)
140         CASE ('U')   ;   prot(:,:) = pyin(:,:) * gcosu(:,:) + pxin(:,:) * gsinu(:,:)
141         CASE ('V')   ;   prot(:,:) = pyin(:,:) * gcosv(:,:) + pxin(:,:) * gsinv(:,:)
142         CASE ('F')   ;   prot(:,:) = pyin(:,:) * gcosf(:,:) + pxin(:,:) * gsinf(:,:)
143         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
144         END SELECT
145      CASE DEFAULT   ;   CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' )
146      END SELECT
147     
148   END SUBROUTINE rot_rep
149
150
151   SUBROUTINE angle
152      !!----------------------------------------------------------------------
153      !!                  ***  ROUTINE angle  ***
154      !!
155      !! ** Purpose :   Compute angles between model grid lines and the North direction
156      !!
157      !! ** Method  :
158      !!
159      !! ** Action  :   Compute (gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf) arrays:
160      !!      sinus and cosinus of the angle between the north-south axe and the
161      !!      j-direction at t, u, v and f-points
162      !!
163      !! History :
164      !!   7.0  !  96-07  (O. Marti )  Original code
165      !!   8.0  !  98-06  (G. Madec )
166      !!   8.5  !  98-06  (G. Madec )  Free form, F90 + opt.
167      !!   9.2  !  07-04  (S. Masson)  Add T, F points and bugfix in cos lateral boundary
168      !!----------------------------------------------------------------------
169      INTEGER ::   ji, jj   ! dummy loop indices
170      INTEGER ::   ierr     ! local integer
171      REAL(wp) ::   &
172         zlam, zphi,            &  ! temporary scalars
173         zlan, zphh,            &  !    "         "
174         zxnpt, zynpt, znnpt,   &  ! x,y components and norm of the vector: T point to North Pole
175         zxnpu, zynpu, znnpu,   &  ! x,y components and norm of the vector: U point to North Pole
176         zxnpv, zynpv, znnpv,   &  ! x,y components and norm of the vector: V point to North Pole
177         zxnpf, zynpf, znnpf,   &  ! x,y components and norm of the vector: F point to North Pole
178         zxvvt, zyvvt, znvvt,   &  ! x,y components and norm of the vector: between V points below and above a T point
179         zxffu, zyffu, znffu,   &  ! x,y components and norm of the vector: between F points below and above a U point
180         zxffv, zyffv, znffv,   &  ! x,y components and norm of the vector: between F points left  and right a V point
181         zxuuf, zyuuf, znuuf       ! x,y components and norm of the vector: between U points below and above a F point
182      !!----------------------------------------------------------------------
183
184      ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj),   & 
185         &      gsinu(jpi,jpj), gcosu(jpi,jpj),   & 
186         &      gsinv(jpi,jpj), gcosv(jpi,jpj),   & 
187         &      gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr )
188      IF(lk_mpp)   CALL mpp_sum( ierr )
189      IF( ierr /= 0 )   CALL ctl_stop('angle: unable to allocate arrays' )
190
191      ! ============================= !
192      ! Compute the cosinus and sinus !
193      ! ============================= !
194      ! (computation done on the north stereographic polar plane)
195
196      DO jj = 2, jpjm1
197!CDIR NOVERRCHK
198         DO ji = fs_2, jpi   ! vector opt.
199
200            ! north pole direction & modulous (at t-point)
201            zlam = glamt(ji,jj)
202            zphi = gphit(ji,jj)
203            zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
204            zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
205            znnpt = zxnpt*zxnpt + zynpt*zynpt
206
207            ! north pole direction & modulous (at u-point)
208            zlam = glamu(ji,jj)
209            zphi = gphiu(ji,jj)
210            zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
211            zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
212            znnpu = zxnpu*zxnpu + zynpu*zynpu
213
214            ! north pole direction & modulous (at v-point)
215            zlam = glamv(ji,jj)
216            zphi = gphiv(ji,jj)
217            zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
218            zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
219            znnpv = zxnpv*zxnpv + zynpv*zynpv
220
221            ! north pole direction & modulous (at f-point)
222            zlam = glamf(ji,jj)
223            zphi = gphif(ji,jj)
224            zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
225            zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
226            znnpf = zxnpf*zxnpf + zynpf*zynpf
227
228            ! j-direction: v-point segment direction (around t-point)
229            zlam = glamv(ji,jj  )
230            zphi = gphiv(ji,jj  )
231            zlan = glamv(ji,jj-1)
232            zphh = gphiv(ji,jj-1)
233            zxvvt =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
234               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
235            zyvvt =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
236               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
237            znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  )
238            znvvt = MAX( znvvt, 1.e-14 )
239
240            ! j-direction: f-point segment direction (around u-point)
241            zlam = glamf(ji,jj  )
242            zphi = gphif(ji,jj  )
243            zlan = glamf(ji,jj-1)
244            zphh = gphif(ji,jj-1)
245            zxffu =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
246               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
247            zyffu =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
248               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
249            znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu )  )
250            znffu = MAX( znffu, 1.e-14 )
251
252            ! i-direction: f-point segment direction (around v-point)
253            zlam = glamf(ji  ,jj)
254            zphi = gphif(ji  ,jj)
255            zlan = glamf(ji-1,jj)
256            zphh = gphif(ji-1,jj)
257            zxffv =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
258               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
259            zyffv =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
260               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
261            znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv )  )
262            znffv = MAX( znffv, 1.e-14 )
263
264            ! j-direction: u-point segment direction (around f-point)
265            zlam = glamu(ji,jj+1)
266            zphi = gphiu(ji,jj+1)
267            zlan = glamu(ji,jj  )
268            zphh = gphiu(ji,jj  )
269            zxuuf =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
270               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
271            zyuuf =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
272               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
273            znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  )
274            znuuf = MAX( znuuf, 1.e-14 )
275
276            ! cosinus and sinus using scalar and vectorial products
277            gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt
278            gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt
279
280            gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu
281            gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu
282
283            gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf
284            gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf
285
286            ! (caution, rotation of 90 degres)
287            gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv
288            gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv
289
290         END DO
291      END DO
292
293      ! =============== !
294      ! Geographic mesh !
295      ! =============== !
296
297      DO jj = 2, jpjm1
298         DO ji = fs_2, jpi   ! vector opt.
299            IF( MOD( ABS( glamv(ji,jj) - glamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
300               gsint(ji,jj) = 0.
301               gcost(ji,jj) = 1.
302            ENDIF
303            IF( MOD( ABS( glamf(ji,jj) - glamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
304               gsinu(ji,jj) = 0.
305               gcosu(ji,jj) = 1.
306            ENDIF
307            IF(      ABS( gphif(ji,jj) - gphif(ji-1,jj) )         < 1.e-8 ) THEN
308               gsinv(ji,jj) = 0.
309               gcosv(ji,jj) = 1.
310            ENDIF
311            IF( MOD( ABS( glamu(ji,jj) - glamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN
312               gsinf(ji,jj) = 0.
313               gcosf(ji,jj) = 1.
314            ENDIF
315         END DO
316      END DO
317
318      ! =========================== !
319      ! Lateral boundary conditions !
320      ! =========================== !
321
322      ! lateral boundary cond.: T-, U-, V-, F-pts, sgn
323      CALL lbc_lnk( gcost, 'T', -1. )   ;   CALL lbc_lnk( gsint, 'T', -1. )
324      CALL lbc_lnk( gcosu, 'U', -1. )   ;   CALL lbc_lnk( gsinu, 'U', -1. )
325      CALL lbc_lnk( gcosv, 'V', -1. )   ;   CALL lbc_lnk( gsinv, 'V', -1. )
326      CALL lbc_lnk( gcosf, 'F', -1. )   ;   CALL lbc_lnk( gsinf, 'F', -1. )
327
328   END SUBROUTINE angle
329
330
331   SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid,     &
332                        pte, ptn )
333      !!----------------------------------------------------------------------
334      !!                    ***  ROUTINE geo2oce  ***
335      !!     
336      !! ** Purpose :
337      !!
338      !! ** Method  :   Change wind stress from geocentric to east/north
339      !!
340      !! History :
341      !!        !         (O. Marti)  Original code
342      !!        !  91-03  (G. Madec)
343      !!        !  92-07  (M. Imbard)
344      !!        !  99-11  (M. Imbard) NetCDF format with IOIPSL
345      !!        !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb
346      !!   8.5  !  02-06  (G. Madec)  F90: Free form
347      !!   3.0  !  07-08  (G. Madec)  geo2oce suppress lon/lat agruments
348      !!----------------------------------------------------------------------
349      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::  pxx, pyy, pzz
350      CHARACTER(len=1)            , INTENT(in   ) ::  cgrid
351      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::  pte, ptn
352      !!
353      REAL(wp), PARAMETER :: rpi = 3.141592653e0
354      REAL(wp), PARAMETER :: rad = rpi / 180.e0
355      INTEGER ::   ig     !
356      INTEGER ::   ierr   ! local integer
357      !!----------------------------------------------------------------------
358
359      IF( .NOT. ALLOCATED( gsinlon ) ) THEN
360         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   &
361            &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
362         IF( lk_mpp    )   CALL mpp_sum( ierr )
363         IF( ierr /= 0 )   CALL ctl_stop('geo2oce: unable to allocate arrays' )
364      ENDIF
365
366      SELECT CASE( cgrid)
367         CASE ( 'T' )   
368            ig = 1
369            IF( .NOT. linit(ig) ) THEN
370               gsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
371               gcoslon(:,:,ig) = COS( rad * glamt(:,:) )
372               gsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
373               gcoslat(:,:,ig) = COS( rad * gphit(:,:) )
374               linit(ig) = .TRUE.
375            ENDIF
376         CASE ( 'U' )   
377            ig = 2
378            IF( .NOT. linit(ig) ) THEN
379               gsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
380               gcoslon(:,:,ig) = COS( rad * glamu(:,:) )
381               gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
382               gcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
383               linit(ig) = .TRUE.
384            ENDIF
385         CASE ( 'V' )   
386            ig = 3
387            IF( .NOT. linit(ig) ) THEN
388               gsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
389               gcoslon(:,:,ig) = COS( rad * glamv(:,:) )
390               gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
391               gcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
392               linit(ig) = .TRUE.
393            ENDIF
394         CASE ( 'F' )   
395            ig = 4
396            IF( .NOT. linit(ig) ) THEN
397               gsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
398               gcoslon(:,:,ig) = COS( rad * glamf(:,:) )
399               gsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
400               gcoslat(:,:,ig) = COS( rad * gphif(:,:) )
401               linit(ig) = .TRUE.
402            ENDIF
403         CASE default   
404            WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
405            CALL ctl_stop( ctmp1 )
406      END SELECT
407     
408      pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy
409      ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx    &
410            - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy    &
411            + gcoslat(:,:,ig) * pzz
412!!$   ptv =   gcoslon(:,:,ig) * gcoslat(:,:,ig) * pxx    &
413!!$         + gsinlon(:,:,ig) * gcoslat(:,:,ig) * pyy    &
414!!$         + gsinlat(:,:,ig) * pzz
415      !
416   END SUBROUTINE geo2oce
417
418   SUBROUTINE oce2geo ( pte, ptn, cgrid,     &
419                        pxx , pyy , pzz )
420      !!----------------------------------------------------------------------
421      !!                    ***  ROUTINE oce2geo  ***
422      !!     
423      !! ** Purpose :
424      !!
425      !! ** Method  :   Change vector from east/north to geocentric
426      !!
427      !! History :
428      !!        !         (A. Caubel)  oce2geo - Original code
429      !!----------------------------------------------------------------------
430      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN    ) ::  pte, ptn
431      CHARACTER(len=1)            , INTENT( IN    ) ::  cgrid
432      REAL(wp), DIMENSION(jpi,jpj), INTENT(   OUT ) ::  pxx , pyy , pzz
433      !!
434      REAL(wp), PARAMETER :: rpi = 3.141592653E0
435      REAL(wp), PARAMETER :: rad = rpi / 180.e0
436      INTEGER ::   ig     !
437      INTEGER ::   ierr   ! local integer
438      !!----------------------------------------------------------------------
439
440      IF( .NOT. ALLOCATED( gsinlon ) ) THEN
441         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   &
442            &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
443         IF( lk_mpp    )   CALL mpp_sum( ierr )
444         IF( ierr /= 0 )   CALL ctl_stop('oce2geo: unable to allocate arrays' )
445      ENDIF
446
447      SELECT CASE( cgrid)
448         CASE ( 'T' )   
449            ig = 1
450            IF( .NOT. linit(ig) ) THEN
451               gsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
452               gcoslon(:,:,ig) = COS( rad * glamt(:,:) )
453               gsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
454               gcoslat(:,:,ig) = COS( rad * gphit(:,:) )
455               linit(ig) = .TRUE.
456            ENDIF
457         CASE ( 'U' )   
458            ig = 2
459            IF( .NOT. linit(ig) ) THEN
460               gsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
461               gcoslon(:,:,ig) = COS( rad * glamu(:,:) )
462               gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
463               gcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
464               linit(ig) = .TRUE.
465            ENDIF
466         CASE ( 'V' )   
467            ig = 3
468            IF( .NOT. linit(ig) ) THEN
469               gsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
470               gcoslon(:,:,ig) = COS( rad * glamv(:,:) )
471               gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
472               gcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
473               linit(ig) = .TRUE.
474            ENDIF
475         CASE ( 'F' )   
476            ig = 4
477            IF( .NOT. linit(ig) ) THEN
478               gsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
479               gcoslon(:,:,ig) = COS( rad * glamf(:,:) )
480               gsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
481               gcoslat(:,:,ig) = COS( rad * gphif(:,:) )
482               linit(ig) = .TRUE.
483            ENDIF
484         CASE default   
485            WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
486            CALL ctl_stop( ctmp1 )
487      END SELECT
488
489       pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn 
490       pyy =   gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn
491       pzz =   gcoslat(:,:,ig) * ptn
492
493     
494   END SUBROUTINE oce2geo
495
496
497   SUBROUTINE repere ( px1, py1, px2, py2, kchoix, cd_type )
498      !!----------------------------------------------------------------------
499      !!                 ***  ROUTINE repere  ***
500      !!       
501      !! ** Purpose :   Change vector componantes between a geopgraphic grid
502      !!      and a stretched coordinates grid.
503      !!
504      !! ** Method  :   
505      !!
506      !! ** Action  :
507      !!
508      !! History :
509      !!        !  89-03  (O. Marti)  original code
510      !!        !  92-02  (M. Imbard)
511      !!        !  93-03  (M. Guyon)  symetrical conditions
512      !!        !  98-05  (B. Blanke)
513      !!   8.5  !  02-08  (G. Madec)  F90: Free form
514      !!----------------------------------------------------------------------
515      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   px1, py1   ! two horizontal components to be rotated
516      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   px2, py2   ! the two horizontal components in the model repere
517      INTEGER , INTENT(in   )                     ::   kchoix     ! type of transformation
518      !                                                           ! = 1 change from geographic to model grid.
519      !                                                           ! =-1 change from model to geographic grid
520      CHARACTER(len=1), INTENT(in   ), OPTIONAL   ::   cd_type    ! define the nature of pt2d array grid-points
521      !
522      CHARACTER(len=1) ::   cl_type      ! define the nature of pt2d array grid-points (T point by default)
523      !!----------------------------------------------------------------------
524
525      cl_type = 'T'
526      IF( PRESENT(cd_type) )   cl_type = cd_type
527         !
528      SELECT CASE (kchoix)
529      CASE ( 1)      ! change from geographic to model grid.
530         CALL rot_rep( px1, py1, cl_type, 'en->i', px2 )
531         CALL rot_rep( px1, py1, cl_type, 'en->j', py2 )
532      CASE (-1)      ! change from model to geographic grid
533         CALL rot_rep( px1, py1, cl_type, 'ij->e', px2 )
534         CALL rot_rep( px1, py1, cl_type, 'ij->n', py2 )
535      CASE DEFAULT   ;   CALL ctl_stop( 'repere: Syntax Error in the definition of kchoix (1 OR -1' )
536      END SELECT
537     
538   END SUBROUTINE repere
539
540
541   SUBROUTINE obs_rot ( psinu, pcosu, psinv, pcosv )
542      !!----------------------------------------------------------------------
543      !!                  ***  ROUTINE obs_rot  ***
544      !!
545      !! ** Purpose :   Copy gsinu, gcosu, gsinv and gsinv
546      !!                to input data for rotations of
547      !!                current at observation points
548      !!
549      !! History :
550      !!   9.2  !  09-02  (K. Mogensen)
551      !!----------------------------------------------------------------------
552      REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT )::   psinu, pcosu, psinv, pcosv   ! copy of data
553      !!----------------------------------------------------------------------
554
555      ! Initialization of gsin* and gcos* at first call
556      ! -----------------------------------------------
557
558      IF( lmust_init ) THEN
559         IF(lwp) WRITE(numout,*)
560         IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched'
561         IF(lwp) WRITE(numout,*) ' ~~~~~~~   coordinate transformation'
562
563         CALL angle       ! initialization of the transformation
564         lmust_init = .FALSE.
565
566      ENDIF
567
568      psinu(:,:) = gsinu(:,:)
569      pcosu(:,:) = gcosu(:,:)
570      psinv(:,:) = gsinv(:,:)
571      pcosv(:,:) = gcosv(:,:)
572
573   END SUBROUTINE obs_rot
574
575  !!======================================================================
576END MODULE geo2ocean
Note: See TracBrowser for help on using the repository browser.