source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90

Last change on this file was 11277, checked in by kingr, 19 months ago

Merged Juan's changes for running AMM15 woth wave coupling.
Corrected minor logic error to allow AMM7-uncoupled to reproduce earlier results.
Few line spacing changes to allow merging with OBS br and trunk rvn 5518.

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