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

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 4 years ago

The Dr Hook changes from my perl code.

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