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/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

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