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 NEMO/trunk/src/OCE/SBC – NEMO

source: NEMO/trunk/src/OCE/SBC/geo2ocean.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

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