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 trunk/NEMO/OPA_SRC – NEMO

source: trunk/NEMO/OPA_SRC/geo2ocean.F90 @ 144

Last change on this file since 144 was 144, checked in by opalod, 20 years ago

CL + CT: BUGFIX089: Syntax correction
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.2 KB
Line 
1MODULE geo2ocean
2   !!======================================================================
3   !!                     ***  MODULE  geo2ocean  ***
4   !! Ocean mesh    :  ???
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !!   repcmo      :
9   !!   angle       :
10   !!   geo2oce     :
11   !!   repere      :   old routine suppress it ???
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE dom_oce         ! mesh and scale factors
15   USE phycst          ! physical constants
16   USE in_out_manager  ! I/O manager
17   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
18
19   IMPLICIT NONE
20
21   !! * Accessibility
22   PRIVATE
23   PUBLIC repcmo   ! routine called by ???.F90
24   PUBLIC geo2oce  ! routine called by ???.F90
25   PUBLIC repere   ! routine called by ???.F90
26
27   !! * Module variables
28   REAL(wp), DIMENSION(jpi,jpj) ::   &
29      gsinu , gcosu ,   &  ! matrix element for change grid u (repcmo.F)
30      gsinv , gcosv ,   &  ! matrix element for change grid v (repcmo.F)
31      gsinus, gcosin      ! matrix element for change grid (repere.F)
32
33  !! * Substitutions
34#  include "vectopt_loop_substitute.h90"
35   !!---------------------------------------------------------------------------------
36   !!   OPA 9.0 , LODYC-IPSL  (2003)
37   !!---------------------------------------------------------------------------------
38
39CONTAINS
40
41   SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1,   &
42                       px2 , py2 , kt )
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE repcmo  ***
45      !!
46      !! ** Purpose :   Change vector componantes from a geographic grid to a
47      !!      stretched coordinates grid.
48      !!
49      !! ** Method  :   Initialization of arrays at the first call.
50      !!
51      !! ** Action  : - px2 : first componante (defined at u point)
52      !!              - py2 : second componante (defined at v point)
53      !!
54      !! History :
55      !!   7.0  !  07-96  (O. Marti)  Original code
56      !!   8.5  !  02-08  (G. Madec)  F90: Free form
57      !!----------------------------------------------------------------------
58      !! * Arguments
59      INTEGER,  INTENT( in ) ::   &
60         kt                ! ocean time-step
61      REAL(wp), INTENT( in ), DIMENSION(jpi,jpj) ::   & 
62         pxu1, pyu1,     & ! geographic vector componantes at u-point
63         pxv1, pyv1        ! geographic vector componantes at v-point
64      REAL(wp), INTENT( out ), DIMENSION(jpi,jpj) ::   & 
65         px2,            & ! i-componante (defined at u-point)
66         py2               ! j-componante (defined at v-point)
67      !!----------------------------------------------------------------------
68
69
70      ! Initialization of gsin* and gcos* at first call
71      ! -----------------------------------------------
72
73      IF( kt <= nit000 + 1 ) THEN
74         IF(lwp) WRITE(numout,*)
75         IF(lwp) WRITE(numout,*) 'repcmo : use the geographic to stretched'
76         IF(lwp) WRITE(numout,*) ' ~~~~~   coordinate transformation'
77
78         CALL angle       ! initialization of the transformation
79      ENDIF
80     
81      ! Change from geographic to stretched coordinate
82      ! ----------------------------------------------
83     
84      px2(:,:) = pxu1(:,:) * gcosu(:,:) + pyu1(:,:) * gsinu(:,:)
85      py2(:,:) = pyv1(:,:) * gcosv(:,:) - pxv1(:,:) * gsinv(:,:)   
86     
87   END SUBROUTINE repcmo
88
89
90   SUBROUTINE angle
91      !!----------------------------------------------------------------------
92      !!                  ***  ROUTINE angle  ***
93      !!
94      !! ** Purpose :   Compute angles between model grid lines and the
95      !!      direction of the North
96      !!
97      !! ** Method  :
98      !!
99      !! ** Action  :   Compute (gsinu, gcosu, gsinv, gcosv) arrays: sinus and
100      !!      cosinus of the angle between the north-south axe and the
101      !!      j-direction at u and v-points
102      !!
103      !! History :
104      !!   7.0  !  96-07  (O. Marti)  Original code
105      !!   8.0  !  98-06  (G. Madec)
106      !!   8.5  !  98-06  (G. Madec)  Free form, F90 + opt.
107      !!----------------------------------------------------------------------
108      !! * local declarations
109      INTEGER ::   ji, jj      ! dummy loop indices
110
111      REAL(wp) ::   &
112         zlam, zphi,             &  ! temporary scalars
113         zlan, zphh,             &  !    "         "
114         zxnpu, zxnpv , znnpu,   &  !    "         "
115         zynpu, zynpv , znnpv,   &  !    "         "
116         zxffu, zmnpfu, zxffv,   &  !    "         "
117         zyffu, zmnpfv, zyffv       !    "         "
118      !!----------------------------------------------------------------------
119
120      ! ============================= !
121      ! Compute the cosinus and sinus !
122      ! ============================= !
123      ! (computation done on the north stereographic polar plane)
124
125      DO jj = 2, jpj
126!CDIR NOVERRCHK
127         DO ji = fs_2, jpi   ! vector opt.
128
129            ! north pole direction & modulous (at u-point)
130            zlam = glamu(ji,jj)
131            zphi = gphiu(ji,jj)
132            zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
133            zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
134            znnpu = zxnpu*zxnpu + zynpu*zynpu
135
136            ! north pole direction & modulous (at v-point)
137            zlam = glamv(ji,jj)
138            zphi = gphiv(ji,jj)
139            zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
140            zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
141            znnpv = zxnpv*zxnpv + zynpv*zynpv
142
143            ! j-direction: f-point segment direction (u-point)
144            zlam = glamf(ji,jj  )
145            zphi = gphif(ji,jj  )
146            zlan = glamf(ji,jj-1)
147            zphh = gphif(ji,jj-1)
148            zxffu =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
149               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
150            zyffu =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
151               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
152            zmnpfu = SQRT ( znnpu * ( zxffu*zxffu + zyffu*zyffu )  )
153            zmnpfu = MAX( zmnpfu, 1.e-14 )
154
155            ! i-direction: f-point segment direction (v-point)
156            zlam = glamf(ji  ,jj)
157            zphi = gphif(ji  ,jj)
158            zlan = glamf(ji-1,jj)
159            zphh = gphif(ji-1,jj)
160            zxffv =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
161               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
162            zyffv =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
163               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
164            zmnpfv = SQRT ( znnpv * ( zxffv*zxffv + zyffv*zyffv )  )
165            zmnpfv = MAX( zmnpfv, 1.e-14 )
166
167            ! cosinus and sinus using scalar and vectorial products
168            gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / zmnpfu
169            gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / zmnpfu
170
171            ! cosinus and sinus using scalar and vectorial products
172            ! (caution, rotation of 90 degres)
173            gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / zmnpfv
174            gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / zmnpfv
175
176         END DO
177      END DO
178
179      ! =============== !
180      ! Geographic mesh !
181      ! =============== !
182
183      DO jj = 2, jpj
184         DO ji = fs_2, jpi   ! vector opt.
185            IF( ABS( glamf(ji,jj) - glamf(ji,jj-1) ) < 1.e-8 ) THEN
186               gsinu(ji,jj) = 0.
187               gcosu(ji,jj) = 1.
188            ENDIF
189            IF( ABS( gphif(ji,jj) - gphif(ji-1,jj) ) < 1.e-8 ) THEN
190               gsinv(ji,jj) = 0.
191               gcosv(ji,jj) = 1.
192            ENDIF
193         END DO
194      END DO
195
196      ! =========================== !
197      ! Lateral boundary conditions !
198      ! =========================== !
199
200      ! lateral boundary cond.: U-, V-pts, sgn
201      CALL lbc_lnk ( gsinu, 'U', -1. )   ;   CALL lbc_lnk( gsinv, 'V', -1. )
202      CALL lbc_lnk ( gcosu, 'U', -1. )   ;   CALL lbc_lnk( gcosv, 'V', -1. )
203
204   END SUBROUTINE angle
205
206
207   SUBROUTINE geo2oce ( pxx , pyy , pzz, cgrid,     &
208                        plon, plat, pte, ptn  , ptv )
209      !!----------------------------------------------------------------------
210      !!                    ***  ROUTINE geo2oce  ***
211      !!     
212      !! ** Purpose :
213      !!
214      !! ** Method  :   Change wind stress from geocentric to east/north
215      !!
216      !! History :
217      !!        !         (O. Marti)  Original code
218      !!        !  91-03  (G. Madec)
219      !!        !  92-07  (M. Imbard)
220      !!        !  99-11  (M. Imbard) NetCDF format with IOIPSL
221      !!        !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb
222      !!   8.5  !  02-06  (G. Madec)  F90: Free form
223      !!----------------------------------------------------------------------
224      !! * Local declarations
225      REAL(wp), INTENT( in ), DIMENSION(jpi,jpj) ::   &
226         pxx, pyy, pzz
227      CHARACTER (len=1), INTENT( in) ::   &
228         cgrid
229      REAL(wp), INTENT( in ), DIMENSION(jpi,jpj) ::   &
230         plon, plat
231      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::    &
232         pte, ptn, ptv
233      REAL(wp), PARAMETER :: rpi = 3.141592653E0
234      REAL(wp), PARAMETER :: rad = rpi / 180.e0
235
236      !! * Local variables
237      INTEGER ::   ig     !
238
239      !! * Local save
240      REAL(wp), SAVE, DIMENSION(jpi,jpj,4) ::   &
241         zsinlon, zcoslon,   &
242         zsinlat, zcoslat
243      LOGICAL, SAVE, DIMENSION (4) ::   &
244         linit = .FALSE.
245      !!----------------------------------------------------------------------
246
247      SELECT CASE( cgrid)
248
249         CASE ( 't' ) ;; ig = 1
250         CASE ( 'u' ) ;; ig = 2
251         CASE ( 'v' ) ;; ig = 3
252         CASE ( 'f' ) ;; ig = 4
253
254         CASE default
255            IF(lwp) WRITE(numout,cform_err)
256            IF(lwp) WRITE(numout,*) 'geo2oce : bad grid argument : ', cgrid
257            nstop = nstop + 1
258       END SELECT
259     
260      IF( .NOT. linit(ig) ) THEN
261         zsinlon (:,:,ig) = SIN (rad * plon)
262         zcoslon (:,:,ig) = COS (rad * plon)
263         zsinlat (:,:,ig) = SIN (rad * plat)
264         zcoslat (:,:,ig) = COS (rad * plat)
265         linit (ig) = .TRUE.
266      ENDIF
267     
268      pte = - zsinlon (:,:,ig) * pxx + zcoslon (:,:,ig) * pyy
269      ptn = - zcoslon (:,:,ig) * zsinlat (:,:,ig) * pxx    &
270            - zsinlon (:,:,ig) * zsinlat (:,:,ig) * pyy    &
271            + zcoslat (:,:,ig) * pzz
272      ptv =   zcoslon (:,:,ig) * zcoslat (:,:,ig) * pxx    &
273            + zsinlon (:,:,ig) * zcoslat (:,:,ig) * pyy    &
274            + zsinlat (:,:,ig) * pzz
275
276   END SUBROUTINE geo2oce
277
278
279   SUBROUTINE repere ( px1, py1, px2, py2, kchoix )
280      !!----------------------------------------------------------------------
281      !!                 ***  ROUTINE repere  ***
282      !!       
283      !! ** Purpose :   Change vector componantes between a geopgraphic grid
284      !!      and a stretched coordinates grid.
285      !!
286      !! ** Method  :   initialization of arrays at the first call.
287      !!
288      !! ** Action  :
289      !!
290      !! History :
291      !!        !  89-03  (O. Marti)  original code
292      !!        !  92-02  (M. Imbard)
293      !!        !  93-03  (M. Guyon)  symetrical conditions
294      !!        !  98-05  (B. Blanke)
295      !!   8.5  !  02-08  (G. Madec)  F90: Free form
296      !!----------------------------------------------------------------------
297      !! * Arguments
298      REAL(wp), INTENT( in   ), DIMENSION(jpi,jpj) ::   &
299         px1, py1          ! two horizontal components to be rotated
300      REAL(wp), INTENT( out  ), DIMENSION(jpi,jpj) ::   &
301         px2, py2          ! the two horizontal components in the model repere
302      INTEGER, INTENT( inout ) ::   &
303         kchoix   ! type of transformation
304                  ! = 1 change from geographic to model grid.
305                  ! =-1 change from model to geographic grid
306                  ! = 0 same as the previous call
307      !! * Local declarations
308      INTEGER, SAVE :: nmem
309
310      INTEGER ::   ji, jj                    ! dummy loop indices
311
312      REAL(wp) :: zxx, zcof1, zcof2,    &
313         ze1t, ze2t
314      REAL(wp), DIMENSION(jpi,jpj) ::   &
315         zlamdu, zphiu,   &
316         zlamdv, zphiv
317      !!----------------------------------------------------------------------
318
319
320      ! 0. Initialization of gsinus and gcosin IF first call
321      ! ----------------------------------------------------
322     
323      ! 0.1 control argument
324     
325      IF( kchoix == 0 ) THEN
326         IF( nmem == 0 ) THEN
327            IF(lwp) WRITE(numout,cform_err)
328            IF(lwp) WRITE(numout,*) 'repere : e r r o r  in kchoix : ', kchoix
329            IF(lwp) WRITE(numout,*) ' for the first call , you must indicate '
330            IF(lwp) WRITE(numout,*) ' the direction of change '
331            IF(lwp) WRITE(numout,*) ' kchoix = 1 geo       --> stretched '
332            IF(lwp) WRITE(numout,*) ' kchoix =-1 stretched --> geo '
333            nstop = nstop + 1
334         ELSE
335            kchoix = nmem
336         ENDIF
337      ELSEIF( kchoix == 1 .OR. kchoix == -1 ) THEN
338         nmem = kchoix
339      ELSE
340         IF(lwp) WRITE(numout,cform_err)
341         IF(lwp) WRITE(numout,*) 'repere : e r r o r  in kchoix : ', kchoix
342         IF(lwp) WRITE(numout,*) ' kchoix must be equal to -1, 0 or 1 '
343         nstop = nstop + 1
344      ENDIF
345
346      ! 0.2 Initialization
347
348      zxx = gsinus(jpi/2,jpj/2)**2+gcosin(jpi/2,jpj/2)**2
349      IF( zxx <= 0.1 ) THEN
350         IF(lwp) WRITE(numout,*) 'repere : initialization '
351         DO jj = 1, jpj
352            DO ji = 2, jpi
353               zlamdu(ji,jj) = glamu(ji,jj) - glamu(ji-1,jj)
354               zlamdu(ji,jj) = ASIN( SIN( rad*zlamdu(ji,jj) ) )/rad
355               zphiu (ji,jj) = gphiu(ji,jj) - gphiu(ji-1,jj)
356            END DO
357         END DO
358         DO jj = 2, jpj
359            zlamdv(:,jj) = glamv(:,jj)-glamv(:,jj-1)
360            zlamdv(:,jj) = ASIN(SIN(rad*zlamdv(:,jj)))/rad
361            zphiv (:,jj) = gphiv(:,jj)-gphiv(:,jj-1)
362         END DO
363         
364         ! 0.3 Boudary conditions and periodicity
365         
366         IF( nperio == 1 .OR.nperio == 4.OR.nperio == 6 ) THEN
367            DO jj = 1, jpj
368               zlamdu(1,jj) = zlamdu(jpim1,jj)
369               zphiu (1,jj) = zphiu (jpim1,jj)
370            END DO
371         ELSE
372            DO jj = 1, jpj
373               zlamdu(1,jj) = zlamdu(2,jj)
374               zphiu (1,jj) = zphiu (2,jj)
375            END DO
376         ENDIF
377         
378         DO ji = 1, jpi
379            zlamdv(ji,1) = zlamdv(ji,2)
380            zphiv (ji,1) = zphiv (ji,2)
381         END DO
382         
383         IF( nperio == 2 ) THEN
384            DO ji = 1, jpi
385               zphiv (ji,1) = zphiv (ji,3)
386            END DO
387         ENDIF
388         
389         ! 0.4 gsinus gcosin
390         
391!CDIR NOVERRCHK
392         DO jj = 1, jpj
393!CDIR NOVERRCHK
394            DO ji = 1, jpi
395               zcof1 = rad * ra * COS( rad * gphit(ji,jj) )
396               zcof2 = rad * ra
397               zlamdu(ji,jj) = zlamdu(ji,jj) * zcof1
398               zlamdv(ji,jj) = zlamdv(ji,jj) * zcof1
399               zphiu (ji,jj) = zphiu (ji,jj) * zcof2
400               zphiv (ji,jj) = zphiv (ji,jj) * zcof2
401            END DO
402         END DO
403
404!CDIR NOVERRCHK
405         DO jj = 1, jpj
406!CDIR NOVERRCHK
407            DO ji = 1, jpi
408               ze1t = SQRT( zlamdu(ji,jj)*zlamdu(ji,jj) + zphiu(ji,jj)*zphiu(ji,jj) )
409               ze2t = SQRT( zlamdv(ji,jj)*zlamdv(ji,jj) + zphiv(ji,jj)*zphiv(ji,jj) )
410               gsinus(ji,jj) = 0.5*( zphiu(ji,jj)/ze1t - zlamdv(ji,jj)/ze2t )
411               gcosin(ji,jj) = 0.5*( zphiv(ji,jj)/ze2t + zlamdu(ji,jj)/ze1t )
412            END DO
413         END DO
414         CALL lbc_lnk( gsinus, 'T', -1. )
415         CALL lbc_lnk( gcosin, 'T', -1. )
416         
417      ENDIF
418
419
420      ! 1. Change from geographic to tretched
421      ! -------------------------------------
422     
423      IF( kchoix == 1 ) THEN
424          px2(:,:) =  px1(:,:) * gcosin(:,:) + py1(:,:) * gsinus(:,:)
425          py2(:,:) = -px1(:,:) * gsinus(:,:) + py1(:,:) * gcosin(:,:)
426      ENDIF
427     
428
429      ! 2. Change from tretched to geographic
430      ! -------------------------------------
431     
432      IF( kchoix == -1 ) THEN
433          px2(:,:) =  px1(:,:) * gcosin(:,:) - py1(:,:) * gsinus(:,:)
434          py2(:,:) =  px1(:,:) * gsinus(:,:) + py1(:,:) * gcosin(:,:)
435      ENDIF
436     
437   END SUBROUTINE repere
438
439  !!======================================================================
440END MODULE geo2ocean
Note: See TracBrowser for help on using the repository browser.