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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 @ 2590

Last change on this file since 2590 was 2590, checked in by trackstand2, 13 years ago

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