source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/geo2ocean.F90 @ 10170

Last change on this file since 10170 was 10170, checked in by smasson, 2 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 2a: add report calls of lbc_lnk, see #2133

  • Property svn:keywords set to Id
File size: 22.1 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  !  06-2006  (G. Madec )  Free form, F90 + opt.
8   !!                 !  04-2007  (S. Masson)  angle: Add T, F points and bugfix in cos lateral boundary
9   !!            3.0  !  07-2008  (G. Madec)  geo2oce suppress lon/lat agruments
10   !!            3.7  !  11-2015  (G. Madec)  remove the unused repere and repcmo routines
11   !!----------------------------------------------------------------------
12#if defined key_agrif
13!clem: these lines do not seem necessary anymore
14!!DIR$ OPTIMIZE (-O 1)  ! cray formulation
15# if defined __INTEL_COMPILER
16!acc: still breaks on at least one Ivybridge cluster with ifort 17.0.4 without this directive
17!DIR$ OPTIMIZE:1        ! intel formulation
18# endif
19#endif
20   !!----------------------------------------------------------------------
21   !!   rot_rep       : Rotate the Repere: geographic grid <==> stretched coordinates grid
22   !!   angle         :
23   !!   geo2oce       :
24   !!   oce2geo       :
25   !!----------------------------------------------------------------------
26   USE dom_oce        ! mesh and scale factors
27   USE phycst         ! physical constants
28   !
29   USE in_out_manager ! I/O manager
30   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
31   USE lib_mpp        ! MPP library
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   rot_rep   ! called in sbccpl, fldread, and cyclone
37   PUBLIC   geo2oce   ! called in sbccpl
38   PUBLIC   oce2geo   ! called in sbccpl
39   PUBLIC   obs_rot   ! called in obs_rot_vel and obs_write
40
41   !                                         ! cos/sin between model grid lines and NP direction
42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsint, gcost   ! at T point
43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinu, gcosu   ! at U point
44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinv, gcosv   ! at V point
45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinf, gcosf   ! at F point
46
47   LOGICAL ,              SAVE, DIMENSION(4)     ::   linit = .FALSE.
48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gsinlon, gcoslon, gsinlat, gcoslat
49
50   LOGICAL ::   lmust_init = .TRUE.        !: used to initialize the cos/sin variables (see above)
51
52   !! * Substitutions
53#  include "vectopt_loop_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
56   !! $Id$
57   !! Software governed by the CeCILL license (see ./LICENSE)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot )
62      !!----------------------------------------------------------------------
63      !!                  ***  ROUTINE rot_rep  ***
64      !!
65      !! ** Purpose :   Rotate the Repere: Change vector componantes between
66      !!                geographic grid <--> stretched coordinates grid.
67      !!----------------------------------------------------------------------
68      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pxin, pyin   ! vector componantes
69      CHARACTER(len=1),             INTENT(in   ) ::   cd_type      ! define the nature of pt2d array grid-points
70      CHARACTER(len=5),             INTENT(in   ) ::   cdtodo       ! type of transpormation:
71      !                                                             ! 'en->i' = east-north to i-component
72      !                                                             ! 'en->j' = east-north to j-component
73      !                                                             ! 'ij->e' = (i,j) components to east
74      !                                                             ! 'ij->n' = (i,j) components to north
75      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   prot     
76      !!----------------------------------------------------------------------
77      !
78      IF( lmust_init ) THEN      ! at 1st call only: set  gsin. & gcos.
79         IF(lwp) WRITE(numout,*)
80         IF(lwp) WRITE(numout,*) ' rot_rep: coordinate transformation : geographic <==> model (i,j)-components'
81         IF(lwp) WRITE(numout,*) ' ~~~~~~~~    '
82         !
83         CALL angle       ! initialization of the transformation
84         lmust_init = .FALSE.
85      ENDIF
86      !
87      SELECT CASE( cdtodo )      ! type of rotation
88      !
89      CASE( 'en->i' )                  ! east-north to i-component
90         SELECT CASE (cd_type)
91         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:)
92         CASE ('U')   ;   prot(:,:) = pxin(:,:) * gcosu(:,:) + pyin(:,:) * gsinu(:,:)
93         CASE ('V')   ;   prot(:,:) = pxin(:,:) * gcosv(:,:) + pyin(:,:) * gsinv(:,:)
94         CASE ('F')   ;   prot(:,:) = pxin(:,:) * gcosf(:,:) + pyin(:,:) * gsinf(:,:)
95         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
96         END SELECT
97      CASE ('en->j')                   ! east-north to j-component
98         SELECT CASE (cd_type)
99         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:)
100         CASE ('U')   ;   prot(:,:) = pyin(:,:) * gcosu(:,:) - pxin(:,:) * gsinu(:,:)
101         CASE ('V')   ;   prot(:,:) = pyin(:,:) * gcosv(:,:) - pxin(:,:) * gsinv(:,:)   
102         CASE ('F')   ;   prot(:,:) = pyin(:,:) * gcosf(:,:) - pxin(:,:) * gsinf(:,:)   
103         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
104         END SELECT
105      CASE ('ij->e')                   ! (i,j)-components to east
106         SELECT CASE (cd_type)
107         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:)
108         CASE ('U')   ;   prot(:,:) = pxin(:,:) * gcosu(:,:) - pyin(:,:) * gsinu(:,:)
109         CASE ('V')   ;   prot(:,:) = pxin(:,:) * gcosv(:,:) - pyin(:,:) * gsinv(:,:)
110         CASE ('F')   ;   prot(:,:) = pxin(:,:) * gcosf(:,:) - pyin(:,:) * gsinf(:,:)
111         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
112         END SELECT
113      CASE ('ij->n')                   ! (i,j)-components to north
114         SELECT CASE (cd_type)
115         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:)
116         CASE ('U')   ;   prot(:,:) = pyin(:,:) * gcosu(:,:) + pxin(:,:) * gsinu(:,:)
117         CASE ('V')   ;   prot(:,:) = pyin(:,:) * gcosv(:,:) + pxin(:,:) * gsinv(:,:)
118         CASE ('F')   ;   prot(:,:) = pyin(:,:) * gcosf(:,:) + pxin(:,:) * gsinf(:,:)
119         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
120         END SELECT
121      CASE DEFAULT   ;   CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' )
122      !
123      END SELECT
124      !
125   END SUBROUTINE rot_rep
126
127
128   SUBROUTINE angle
129      !!----------------------------------------------------------------------
130      !!                  ***  ROUTINE angle  ***
131      !!
132      !! ** Purpose :   Compute angles between model grid lines and the North direction
133      !!
134      !! ** Method  :   sinus and cosinus of the angle between the north-south axe
135      !!              and the j-direction at t, u, v and f-points
136      !!                dot and cross products are used to obtain cos and sin, resp.
137      !!
138      !! ** Action  : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf
139      !!----------------------------------------------------------------------
140      INTEGER  ::   ji, jj   ! dummy loop indices
141      INTEGER  ::   ierr     ! local integer
142      REAL(wp) ::   zlam, zphi            ! local scalars
143      REAL(wp) ::   zlan, zphh            !   -      -
144      REAL(wp) ::   zxnpt, zynpt, znnpt   ! x,y components and norm of the vector: T point to North Pole
145      REAL(wp) ::   zxnpu, zynpu, znnpu   ! x,y components and norm of the vector: U point to North Pole
146      REAL(wp) ::   zxnpv, zynpv, znnpv   ! x,y components and norm of the vector: V point to North Pole
147      REAL(wp) ::   zxnpf, zynpf, znnpf   ! x,y components and norm of the vector: F point to North Pole
148      REAL(wp) ::   zxvvt, zyvvt, znvvt   ! x,y components and norm of the vector: between V points below and above a T point
149      REAL(wp) ::   zxffu, zyffu, znffu   ! x,y components and norm of the vector: between F points below and above a U point
150      REAL(wp) ::   zxffv, zyffv, znffv   ! x,y components and norm of the vector: between F points left  and right a V point
151      REAL(wp) ::   zxuuf, zyuuf, znuuf   ! x,y components and norm of the vector: between U points below and above a F point
152      !!----------------------------------------------------------------------
153      !
154      ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj),   & 
155         &      gsinu(jpi,jpj), gcosu(jpi,jpj),   & 
156         &      gsinv(jpi,jpj), gcosv(jpi,jpj),   & 
157         &      gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr )
158      IF(lk_mpp)   CALL mpp_sum( ierr )
159      IF( ierr /= 0 )   CALL ctl_stop( 'angle: unable to allocate arrays' )
160      !
161      ! ============================= !
162      ! Compute the cosinus and sinus !
163      ! ============================= !
164      ! (computation done on the north stereographic polar plane)
165      !
166      DO jj = 2, jpjm1
167         DO ji = fs_2, jpi   ! vector opt.
168            !                 
169            zlam = glamt(ji,jj)     ! north pole direction & modulous (at t-point)
170            zphi = gphit(ji,jj)
171            zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
172            zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
173            znnpt = zxnpt*zxnpt + zynpt*zynpt
174            !
175            zlam = glamu(ji,jj)     ! north pole direction & modulous (at u-point)
176            zphi = gphiu(ji,jj)
177            zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
178            zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
179            znnpu = zxnpu*zxnpu + zynpu*zynpu
180            !
181            zlam = glamv(ji,jj)     ! north pole direction & modulous (at v-point)
182            zphi = gphiv(ji,jj)
183            zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
184            zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
185            znnpv = zxnpv*zxnpv + zynpv*zynpv
186            !
187            zlam = glamf(ji,jj)     ! north pole direction & modulous (at f-point)
188            zphi = gphif(ji,jj)
189            zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
190            zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
191            znnpf = zxnpf*zxnpf + zynpf*zynpf
192            !
193            zlam = glamv(ji,jj  )   ! j-direction: v-point segment direction (around t-point)
194            zphi = gphiv(ji,jj  )
195            zlan = glamv(ji,jj-1)
196            zphh = gphiv(ji,jj-1)
197            zxvvt =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
198               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
199            zyvvt =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
200               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
201            znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  )
202            znvvt = MAX( znvvt, 1.e-14 )
203            !
204            zlam = glamf(ji,jj  )   ! j-direction: f-point segment direction (around u-point)
205            zphi = gphif(ji,jj  )
206            zlan = glamf(ji,jj-1)
207            zphh = gphif(ji,jj-1)
208            zxffu =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
209               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
210            zyffu =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
211               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
212            znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu )  )
213            znffu = MAX( znffu, 1.e-14 )
214            !
215            zlam = glamf(ji  ,jj)   ! i-direction: f-point segment direction (around v-point)
216            zphi = gphif(ji  ,jj)
217            zlan = glamf(ji-1,jj)
218            zphh = gphif(ji-1,jj)
219            zxffv =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
220               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
221            zyffv =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
222               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
223            znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv )  )
224            znffv = MAX( znffv, 1.e-14 )
225            !
226            zlam = glamu(ji,jj+1)   ! j-direction: u-point segment direction (around f-point)
227            zphi = gphiu(ji,jj+1)
228            zlan = glamu(ji,jj  )
229            zphh = gphiu(ji,jj  )
230            zxuuf =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
231               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
232            zyuuf =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
233               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
234            znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  )
235            znuuf = MAX( znuuf, 1.e-14 )
236            !
237            !                       ! cosinus and sinus using dot and cross products
238            gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt
239            gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt
240            !
241            gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu
242            gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu
243            !
244            gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf
245            gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf
246            !
247            gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv
248            gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv     ! (caution, rotation of 90 degres)
249            !
250         END DO
251      END DO
252
253      ! =============== !
254      ! Geographic mesh !
255      ! =============== !
256
257      DO jj = 2, jpjm1
258         DO ji = fs_2, jpi   ! vector opt.
259            IF( MOD( ABS( glamv(ji,jj) - glamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
260               gsint(ji,jj) = 0.
261               gcost(ji,jj) = 1.
262            ENDIF
263            IF( MOD( ABS( glamf(ji,jj) - glamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
264               gsinu(ji,jj) = 0.
265               gcosu(ji,jj) = 1.
266            ENDIF
267            IF(      ABS( gphif(ji,jj) - gphif(ji-1,jj) )         < 1.e-8 ) THEN
268               gsinv(ji,jj) = 0.
269               gcosv(ji,jj) = 1.
270            ENDIF
271            IF( MOD( ABS( glamu(ji,jj) - glamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN
272               gsinf(ji,jj) = 0.
273               gcosf(ji,jj) = 1.
274            ENDIF
275         END DO
276      END DO
277
278      ! =========================== !
279      ! Lateral boundary conditions !
280      ! =========================== !
281      !           ! lateral boundary cond.: T-, U-, V-, F-pts, sgn
282      CALL lbc_lnk_multi( 'geo2ocean', gcost, 'T', -1., gsint, 'T', -1., gcosu, 'U', -1., gsinu, 'U', -1., & 
283                      &   gcosv, 'V', -1., gsinv, 'V', -1., gcosf, 'F', -1., gsinf, 'F', -1.  )
284      !
285   END SUBROUTINE angle
286
287
288   SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, pte, ptn )
289      !!----------------------------------------------------------------------
290      !!                    ***  ROUTINE geo2oce  ***
291      !!     
292      !! ** Purpose :
293      !!
294      !! ** Method  :   Change a vector from geocentric to east/north
295      !!
296      !!----------------------------------------------------------------------
297      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::  pxx, pyy, pzz
298      CHARACTER(len=1)            , INTENT(in   ) ::  cgrid
299      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::  pte, ptn
300      !
301      REAL(wp), PARAMETER :: rpi = 3.141592653e0
302      REAL(wp), PARAMETER :: rad = rpi / 180.e0
303      INTEGER ::   ig     !
304      INTEGER ::   ierr   ! local integer
305      !!----------------------------------------------------------------------
306      !
307      IF( .NOT. ALLOCATED( gsinlon ) ) THEN
308         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   &
309            &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
310         IF( lk_mpp    )   CALL mpp_sum( ierr )
311         IF( ierr /= 0 )   CALL ctl_stop('geo2oce: unable to allocate arrays' )
312      ENDIF
313      !
314      SELECT CASE( cgrid)
315      CASE ( 'T' )   
316         ig = 1
317         IF( .NOT. linit(ig) ) THEN
318            gsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
319            gcoslon(:,:,ig) = COS( rad * glamt(:,:) )
320            gsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
321            gcoslat(:,:,ig) = COS( rad * gphit(:,:) )
322            linit(ig) = .TRUE.
323         ENDIF
324      CASE ( 'U' )   
325         ig = 2
326         IF( .NOT. linit(ig) ) THEN
327            gsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
328            gcoslon(:,:,ig) = COS( rad * glamu(:,:) )
329            gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
330            gcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
331            linit(ig) = .TRUE.
332         ENDIF
333      CASE ( 'V' )   
334         ig = 3
335         IF( .NOT. linit(ig) ) THEN
336            gsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
337            gcoslon(:,:,ig) = COS( rad * glamv(:,:) )
338            gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
339            gcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
340            linit(ig) = .TRUE.
341         ENDIF
342      CASE ( 'F' )   
343         ig = 4
344         IF( .NOT. linit(ig) ) THEN
345            gsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
346            gcoslon(:,:,ig) = COS( rad * glamf(:,:) )
347            gsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
348            gcoslat(:,:,ig) = COS( rad * gphif(:,:) )
349            linit(ig) = .TRUE.
350         ENDIF
351      CASE default   
352         WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
353         CALL ctl_stop( ctmp1 )
354      END SELECT
355      !
356      pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy
357      ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx    &
358         &  - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy    &
359         &  + gcoslat(:,:,ig) * pzz
360      !
361   END SUBROUTINE geo2oce
362
363
364   SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz )
365      !!----------------------------------------------------------------------
366      !!                    ***  ROUTINE oce2geo  ***
367      !!     
368      !! ** Purpose :
369      !!
370      !! ** Method  :   Change vector from east/north to geocentric
371      !!
372      !! History :     ! (A. Caubel)  oce2geo - Original code
373      !!----------------------------------------------------------------------
374      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN    ) ::  pte, ptn
375      CHARACTER(len=1)            , INTENT( IN    ) ::  cgrid
376      REAL(wp), DIMENSION(jpi,jpj), INTENT(   OUT ) ::  pxx , pyy , pzz
377      !!
378      REAL(wp), PARAMETER :: rpi = 3.141592653E0
379      REAL(wp), PARAMETER :: rad = rpi / 180.e0
380      INTEGER ::   ig     !
381      INTEGER ::   ierr   ! local integer
382      !!----------------------------------------------------------------------
383
384      IF( .NOT. ALLOCATED( gsinlon ) ) THEN
385         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   &
386            &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
387         IF( lk_mpp    )   CALL mpp_sum( ierr )
388         IF( ierr /= 0 )   CALL ctl_stop('oce2geo: unable to allocate arrays' )
389      ENDIF
390
391      SELECT CASE( cgrid)
392         CASE ( 'T' )   
393            ig = 1
394            IF( .NOT. linit(ig) ) THEN
395               gsinlon(:,:,ig) = SIN( rad * glamt(:,:) )
396               gcoslon(:,:,ig) = COS( rad * glamt(:,:) )
397               gsinlat(:,:,ig) = SIN( rad * gphit(:,:) )
398               gcoslat(:,:,ig) = COS( rad * gphit(:,:) )
399               linit(ig) = .TRUE.
400            ENDIF
401         CASE ( 'U' )   
402            ig = 2
403            IF( .NOT. linit(ig) ) THEN
404               gsinlon(:,:,ig) = SIN( rad * glamu(:,:) )
405               gcoslon(:,:,ig) = COS( rad * glamu(:,:) )
406               gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) )
407               gcoslat(:,:,ig) = COS( rad * gphiu(:,:) )
408               linit(ig) = .TRUE.
409            ENDIF
410         CASE ( 'V' )   
411            ig = 3
412            IF( .NOT. linit(ig) ) THEN
413               gsinlon(:,:,ig) = SIN( rad * glamv(:,:) )
414               gcoslon(:,:,ig) = COS( rad * glamv(:,:) )
415               gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) )
416               gcoslat(:,:,ig) = COS( rad * gphiv(:,:) )
417               linit(ig) = .TRUE.
418            ENDIF
419         CASE ( 'F' )   
420            ig = 4
421            IF( .NOT. linit(ig) ) THEN
422               gsinlon(:,:,ig) = SIN( rad * glamf(:,:) )
423               gcoslon(:,:,ig) = COS( rad * glamf(:,:) )
424               gsinlat(:,:,ig) = SIN( rad * gphif(:,:) )
425               gcoslat(:,:,ig) = COS( rad * gphif(:,:) )
426               linit(ig) = .TRUE.
427            ENDIF
428         CASE default   
429            WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
430            CALL ctl_stop( ctmp1 )
431      END SELECT
432      !
433      pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn 
434      pyy =   gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn
435      pzz =   gcoslat(:,:,ig) * ptn
436      !
437   END SUBROUTINE oce2geo
438
439
440   SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv )
441      !!----------------------------------------------------------------------
442      !!                  ***  ROUTINE obs_rot  ***
443      !!
444      !! ** Purpose :   Copy gsinu, gcosu, gsinv and gsinv
445      !!                to input data for rotations of
446      !!                current at observation points
447      !!
448      !! History :  9.2  !  09-02  (K. Mogensen)
449      !!----------------------------------------------------------------------
450      REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT )::   psinu, pcosu, psinv, pcosv   ! copy of data
451      !!----------------------------------------------------------------------
452      !
453      ! Initialization of gsin* and gcos* at first call
454      ! -----------------------------------------------
455      IF( lmust_init ) THEN
456         IF(lwp) WRITE(numout,*)
457         IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched'
458         IF(lwp) WRITE(numout,*) ' ~~~~~~~   coordinate transformation'
459         CALL angle       ! initialization of the transformation
460         lmust_init = .FALSE.
461      ENDIF
462      !
463      psinu(:,:) = gsinu(:,:)
464      pcosu(:,:) = gcosu(:,:)
465      psinv(:,:) = gsinv(:,:)
466      pcosv(:,:) = gcosv(:,:)
467      !
468   END SUBROUTINE obs_rot
469
470  !!======================================================================
471END MODULE geo2ocean
Note: See TracBrowser for help on using the repository browser.