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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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