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 @ 88

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

CT : UPDATE057 : # General syntax, alignement, comments corrections

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