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/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/geo2ocean.F90 @ 12340

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

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

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