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

source: branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 @ 6251

Last change on this file since 6251 was 6251, checked in by frrh, 8 years ago

Merge fcm:NEMO.xm/branches/UKMO/dev_r5107_hadgem3_cplfld@5592

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