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.
Changeset 672 – NEMO

Changeset 672


Ignore:
Timestamp:
2007-06-01T19:35:08+02:00 (17 years ago)
Author:
ctlod
Message:

nemo_v2_bugfix_043 : SM : bugfix on lbc + add new features

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/geo2ocean.F90

    r474 r672  
    2121   !! * Accessibility 
    2222   PRIVATE 
    23    PUBLIC repcmo   ! routine called by ???.F90 
    24    PUBLIC geo2oce  ! routine called by ???.F90 
    25    PUBLIC repere   ! routine called by ???.F90 
     23   PUBLIC rot_rep, repcmo, repere, geo2oce   ! only rot_rep should be used 
     24                                             ! repcmo and repere are keep only for compatibility. 
     25                                             ! they are only a useless overlay of rot_rep 
    2626 
    2727   !! * Module variables 
    2828   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) 
     29      gsint, gcost,   &  ! cos/sin between model grid lines and NP direction at T point 
     30      gsinu, gcosu,   &  ! cos/sin between model grid lines and NP direction at U point 
     31      gsinv, gcosv,   &  ! cos/sin between model grid lines and NP direction at V point 
     32      gsinf, gcosf       ! cos/sin between model grid lines and NP direction at F point 
     33 
     34   LOGICAL ::   lmust_init = .TRUE.        !: used to initialize the cos/sin variables (se above) 
    3235 
    3336  !! * Substitutions 
     
    6871         py2               ! j-componante (defined at v-point) 
    6972      !!---------------------------------------------------------------------- 
    70  
     73       
     74      ! Change from geographic to stretched coordinate 
     75      ! ---------------------------------------------- 
     76       
     77      px2(:,:) = rot_rep( pxu1, pyu1, 'U', 'en->i' ) 
     78      py2(:,:) = rot_rep( pxv1, pyv1, 'V', 'en->j' ) 
     79       
     80   END SUBROUTINE repcmo 
     81 
     82 
     83   FUNCTION rot_rep ( pxin, pyin, cd_type, cdtodo ) 
     84      !!---------------------------------------------------------------------- 
     85      !!                  ***  ROUTINE rot_rep  *** 
     86      !! 
     87      !! ** Purpose :   Rotate the Repere: Change vector componantes between 
     88      !!                geographic grid <--> stretched coordinates grid. 
     89      !! 
     90      !! History : 
     91      !!   9.2  !  07-04  (S. Masson)   
     92      !!                  (O. Marti ) Original code (repere and repcmo) 
     93      !!---------------------------------------------------------------------- 
     94      !! * Arguments  
     95      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) ::   pxin, pyin   ! vector componantes 
     96      CHARACTER(len=1),             INTENT( IN ) ::   cd_type      ! define the nature of pt2d array grid-points 
     97      CHARACTER(len=5),             INTENT( IN ) ::   cdtodo       ! specify the work to do: 
     98      !!                                                           ! 'en->i' east-north componantes to model i componante 
     99      !!                                                           ! 'en->j' east-north componantes to model j componante 
     100      !!                                                           ! 'ij->e' model i-j componantes to east componante 
     101      !!                                                           ! 'ij->n' model i-j componantes to east componante 
     102      REAL(wp), DIMENSION(jpi,jpj) ::   rot_rep       
     103 
     104      !!---------------------------------------------------------------------- 
    71105 
    72106      ! Initialization of gsin* and gcos* at first call 
    73107      ! ----------------------------------------------- 
    74108 
    75       IF( kt <= nit000 + 1 ) THEN 
     109      IF( lmust_init ) THEN 
    76110         IF(lwp) WRITE(numout,*) 
    77          IF(lwp) WRITE(numout,*) 'repcmo : use the geographic to stretched' 
    78          IF(lwp) WRITE(numout,*) ' ~~~~~   coordinate transformation' 
     111         IF(lwp) WRITE(numout,*) ' rot_rep : geographic <--> stretched' 
     112         IF(lwp) WRITE(numout,*) ' ~~~~~    coordinate transformation' 
    79113 
    80114         CALL angle       ! initialization of the transformation 
     115         lmust_init = .FALSE. 
     116 
    81117      ENDIF 
    82118       
    83       ! Change from geographic to stretched coordinate 
    84       ! ---------------------------------------------- 
    85        
    86       px2(:,:) = pxu1(:,:) * gcosu(:,:) + pyu1(:,:) * gsinu(:,:) 
    87       py2(:,:) = pyv1(:,:) * gcosv(:,:) - pxv1(:,:) * gsinv(:,:)    
    88        
    89    END SUBROUTINE repcmo 
     119      SELECT CASE (cdtodo) 
     120      CASE ('en->i')      ! 'en->i' est-north componantes to model i componante 
     121         SELECT CASE (cd_type) 
     122         CASE ('T')   ;   rot_rep(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:) 
     123         CASE ('U')   ;   rot_rep(:,:) = pxin(:,:) * gcosu(:,:) + pyin(:,:) * gsinu(:,:) 
     124         CASE ('V')   ;   rot_rep(:,:) = pxin(:,:) * gcosv(:,:) + pyin(:,:) * gsinv(:,:) 
     125         CASE ('F')   ;   rot_rep(:,:) = pxin(:,:) * gcosf(:,:) + pyin(:,:) * gsinf(:,:) 
     126         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 
     127         END SELECT 
     128      CASE ('en->j')      ! 'en->j' est-north componantes to model j componante 
     129         SELECT CASE (cd_type) 
     130         CASE ('T')   ;   rot_rep(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:) 
     131         CASE ('U')   ;   rot_rep(:,:) = pyin(:,:) * gcosu(:,:) - pxin(:,:) * gsinu(:,:) 
     132         CASE ('V')   ;   rot_rep(:,:) = pyin(:,:) * gcosv(:,:) - pxin(:,:) * gsinv(:,:)    
     133         CASE ('F')   ;   rot_rep(:,:) = pyin(:,:) * gcosf(:,:) - pxin(:,:) * gsinf(:,:)    
     134         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 
     135         END SELECT 
     136      CASE ('ij->e')      ! 'ij->e' model i-j componantes to est componante 
     137         SELECT CASE (cd_type) 
     138         CASE ('T')   ;   rot_rep(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:) 
     139         CASE ('U')   ;   rot_rep(:,:) = pxin(:,:) * gcosu(:,:) - pyin(:,:) * gsinu(:,:) 
     140         CASE ('V')   ;   rot_rep(:,:) = pxin(:,:) * gcosv(:,:) - pyin(:,:) * gsinv(:,:) 
     141         CASE ('F')   ;   rot_rep(:,:) = pxin(:,:) * gcosf(:,:) - pyin(:,:) * gsinf(:,:) 
     142         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 
     143         END SELECT 
     144      CASE ('ij->n')      ! 'ij->n' model i-j componantes to est componante 
     145         SELECT CASE (cd_type) 
     146         CASE ('T')   ;   rot_rep(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:) 
     147         CASE ('U')   ;   rot_rep(:,:) = pyin(:,:) * gcosu(:,:) + pxin(:,:) * gsinu(:,:) 
     148         CASE ('V')   ;   rot_rep(:,:) = pyin(:,:) * gcosv(:,:) + pxin(:,:) * gsinv(:,:) 
     149         CASE ('F')   ;   rot_rep(:,:) = pyin(:,:) * gcosf(:,:) + pxin(:,:) * gsinf(:,:) 
     150         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 
     151         END SELECT 
     152      CASE DEFAULT   ;   CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' ) 
     153      END SELECT 
     154       
     155   END FUNCTION rot_rep 
    90156 
    91157 
     
    94160      !!                  ***  ROUTINE angle  *** 
    95161      !!  
    96       !! ** Purpose :   Compute angles between model grid lines and the  
    97       !!      direction of the North 
     162      !! ** Purpose :   Compute angles between model grid lines and the North direction 
    98163      !! 
    99164      !! ** Method  : 
    100165      !! 
    101       !! ** Action  :   Compute (gsinu, gcosu, gsinv, gcosv) arrays: sinus and  
    102       !!      cosinus of the angle between the north-south axe and the  
    103       !!      j-direction at u and v-points 
     166      !! ** Action  :   Compute (gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf) arrays: 
     167      !!      sinus and cosinus of the angle between the north-south axe and the  
     168      !!      j-direction at t, u, v and f-points 
    104169      !! 
    105170      !! History : 
    106       !!   7.0  !  96-07  (O. Marti)  Original code 
    107       !!   8.0  !  98-06  (G. Madec) 
    108       !!   8.5  !  98-06  (G. Madec)  Free form, F90 + opt. 
     171      !!   7.0  !  96-07  (O. Marti )  Original code 
     172      !!   8.0  !  98-06  (G. Madec ) 
     173      !!   8.5  !  98-06  (G. Madec )  Free form, F90 + opt. 
     174      !!   9.2  !  07-04  (S. Masson)  Add T, F points and bugfix in cos lateral boundary 
    109175      !!---------------------------------------------------------------------- 
    110176      !! * local declarations 
     
    112178 
    113179      REAL(wp) ::   & 
    114          zlam, zphi,             &  ! temporary scalars 
    115          zlan, zphh,             &  !    "         " 
    116          zxnpu, zxnpv , znnpu,   &  !    "         " 
    117          zynpu, zynpv , znnpv,   &  !    "         " 
    118          zxffu, zmnpfu, zxffv,   &  !    "         " 
    119          zyffu, zmnpfv, zyffv       !    "         " 
     180         zlam, zphi,            &  ! temporary scalars 
     181         zlan, zphh,            &  !    "         " 
     182         zxnpt, zynpt, znnpt,   &  ! x,y components and norm of the vector: T point to North Pole 
     183         zxnpu, zynpu, znnpu,   &  ! x,y components and norm of the vector: U point to North Pole 
     184         zxnpv, zynpv, znnpv,   &  ! x,y components and norm of the vector: V point to North Pole 
     185         zxnpf, zynpf, znnpf,   &  ! x,y components and norm of the vector: F point to North Pole 
     186         zxvvt, zyvvt, znvvt,   &  ! x,y components and norm of the vector: between V points below and above a T point 
     187         zxffu, zyffu, znffu,   &  ! x,y components and norm of the vector: between F points below and above a U point 
     188         zxffv, zyffv, znffv,   &  ! x,y components and norm of the vector: between F points left  and right a V point 
     189         zxuuf, zyuuf, znuuf       ! x,y components and norm of the vector: between U points below and above a F point 
    120190      !!---------------------------------------------------------------------- 
    121191 
     
    125195      ! (computation done on the north stereographic polar plane) 
    126196 
    127       DO jj = 2, jpj 
     197      DO jj = 2, jpjm1 
    128198!CDIR NOVERRCHK 
    129199         DO ji = fs_2, jpi   ! vector opt. 
     200 
     201            ! north pole direction & modulous (at t-point) 
     202            zlam = glamt(ji,jj) 
     203            zphi = gphit(ji,jj) 
     204            zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
     205            zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
     206            znnpt = zxnpt*zxnpt + zynpt*zynpt 
    130207 
    131208            ! north pole direction & modulous (at u-point) 
     
    143220            znnpv = zxnpv*zxnpv + zynpv*zynpv 
    144221 
    145             ! j-direction: f-point segment direction (u-point) 
     222            ! north pole direction & modulous (at f-point) 
     223            zlam = glamf(ji,jj) 
     224            zphi = gphif(ji,jj) 
     225            zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
     226            zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
     227            znnpf = zxnpf*zxnpf + zynpf*zynpf 
     228 
     229            ! j-direction: v-point segment direction (around t-point) 
     230            zlam = glamv(ji,jj  ) 
     231            zphi = gphiv(ji,jj  ) 
     232            zlan = glamv(ji,jj-1) 
     233            zphh = gphiv(ji,jj-1) 
     234            zxvvt =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   & 
     235               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 
     236            zyvvt =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   & 
     237               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 
     238            znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  ) 
     239            znvvt = MAX( znvvt, 1.e-14 ) 
     240 
     241            ! j-direction: f-point segment direction (around u-point) 
    146242            zlam = glamf(ji,jj  ) 
    147243            zphi = gphif(ji,jj  ) 
     
    152248            zyffu =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   & 
    153249               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 
    154             zmnpfu = SQRT ( znnpu * ( zxffu*zxffu + zyffu*zyffu )  ) 
    155             zmnpfu = MAX( zmnpfu, 1.e-14 ) 
    156  
    157             ! i-direction: f-point segment direction (v-point) 
     250            znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu )  ) 
     251            znffu = MAX( znffu, 1.e-14 ) 
     252 
     253            ! i-direction: f-point segment direction (around v-point) 
    158254            zlam = glamf(ji  ,jj) 
    159255            zphi = gphif(ji  ,jj) 
     
    164260            zyffv =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   & 
    165261               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 
    166             zmnpfv = SQRT ( znnpv * ( zxffv*zxffv + zyffv*zyffv )  ) 
    167             zmnpfv = MAX( zmnpfv, 1.e-14 ) 
     262            znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv )  ) 
     263            znffv = MAX( znffv, 1.e-14 ) 
     264 
     265            ! j-direction: u-point segment direction (around f-point) 
     266            zlam = glamu(ji,jj+1) 
     267            zphi = gphiu(ji,jj+1) 
     268            zlan = glamu(ji,jj  ) 
     269            zphh = gphiu(ji,jj  ) 
     270            zxuuf =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   & 
     271               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 
     272            zyuuf =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   & 
     273               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 
     274            znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  ) 
     275            znuuf = MAX( znuuf, 1.e-14 ) 
    168276 
    169277            ! cosinus and sinus using scalar and vectorial products 
    170             gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / zmnpfu 
    171             gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / zmnpfu 
    172  
    173             ! cosinus and sinus using scalar and vectorial products 
     278            gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt 
     279            gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt 
     280 
     281            gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu 
     282            gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu 
     283 
     284            gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf 
     285            gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf 
     286 
    174287            ! (caution, rotation of 90 degres) 
    175             gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / zmnpfv 
    176             gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / zmnpfv 
     288            gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv 
     289            gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv 
    177290 
    178291         END DO 
     
    183296      ! =============== ! 
    184297 
    185       DO jj = 2, jpj 
     298      DO jj = 2, jpjm1 
    186299         DO ji = fs_2, jpi   ! vector opt. 
    187             IF( ABS( glamf(ji,jj) - glamf(ji,jj-1) ) < 1.e-8 ) THEN 
     300            IF( MOD( ABS( glamv(ji,jj) - glamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN 
     301               gsint(ji,jj) = 0. 
     302               gcost(ji,jj) = 1. 
     303            ENDIF 
     304            IF( MOD( ABS( glamf(ji,jj) - glamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN 
    188305               gsinu(ji,jj) = 0. 
    189306               gcosu(ji,jj) = 1. 
    190307            ENDIF 
    191             IF( ABS( gphif(ji,jj) - gphif(ji-1,jj) ) < 1.e-8 ) THEN 
     308            IF(      ABS( gphif(ji,jj) - gphif(ji-1,jj) )        < 1.e-8 ) THEN 
    192309               gsinv(ji,jj) = 0. 
    193310               gcosv(ji,jj) = 1. 
    194311            ENDIF 
     312            IF( MOD( ABS( glamu(ji,jj) - glamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN 
     313               gsinf(ji,jj) = 0. 
     314               gcosf(ji,jj) = 1. 
     315            ENDIF 
    195316         END DO 
    196317      END DO 
     
    200321      ! =========================== ! 
    201322 
    202       ! lateral boundary cond.: U-, V-pts, sgn 
    203       CALL lbc_lnk ( gsinu, 'U', -1. )   ;   CALL lbc_lnk( gsinv, 'V', -1. ) 
    204       CALL lbc_lnk ( gcosu, 'U', -1. )   ;   CALL lbc_lnk( gcosv, 'V', -1. ) 
     323      ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 
     324      CALL lbc_lnk ( gcost, 'T', 1. )   ;   CALL lbc_lnk( gsint, 'T', -1. ) 
     325      CALL lbc_lnk ( gcosu, 'U', 1. )   ;   CALL lbc_lnk( gsinu, 'U', -1. ) 
     326      CALL lbc_lnk ( gcosv, 'V', 1. )   ;   CALL lbc_lnk( gsinv, 'V', -1. ) 
     327      CALL lbc_lnk ( gcosf, 'F', 1. )   ;   CALL lbc_lnk( gsinf, 'F', -1. ) 
    205328 
    206329   END SUBROUTINE angle 
     
    278401 
    279402 
    280    SUBROUTINE repere ( px1, py1, px2, py2, kchoix ) 
     403   SUBROUTINE repere ( px1, py1, px2, py2, kchoix, cd_type ) 
    281404      !!---------------------------------------------------------------------- 
    282405      !!                 ***  ROUTINE repere  *** 
     
    285408      !!      and a stretched coordinates grid. 
    286409      !! 
    287       !! ** Method  :   initialization of arrays at the first call. 
     410      !! ** Method  :    
    288411      !! 
    289412      !! ** Action  : 
     
    297420      !!---------------------------------------------------------------------- 
    298421      !! * Arguments 
    299       REAL(wp), INTENT( in   ), DIMENSION(jpi,jpj) ::   & 
     422      REAL(wp), INTENT( IN   ), DIMENSION(jpi,jpj) ::   & 
    300423         px1, py1          ! two horizontal components to be rotated 
    301       REAL(wp), INTENT( out  ), DIMENSION(jpi,jpj) ::   & 
     424      REAL(wp), INTENT( OUT  ), DIMENSION(jpi,jpj) ::   & 
    302425         px2, py2          ! the two horizontal components in the model repere 
    303       INTEGER, INTENT( inout ) ::   & 
     426      INTEGER, INTENT( IN ) ::   & 
    304427         kchoix   ! type of transformation 
    305428                  ! = 1 change from geographic to model grid. 
    306429                  ! =-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             WRITE(ctmp1,*) 'repere : e r r o r  in kchoix : ', kchoix 
    329             CALL ctl_stop( ctmp1, ' for the first call , you must indicate ',   & 
    330                  &                ' the direction of change ',   & 
    331                  &                ' kchoix = 1 geo       --> stretched ',   & 
    332                  &                ' kchoix =-1 stretched --> geo ' ) 
    333          ELSE 
    334             kchoix = nmem 
    335          ENDIF 
    336       ELSEIF( kchoix == 1 .OR. kchoix == -1 ) THEN 
    337          nmem = kchoix 
    338       ELSE 
    339          WRITE(ctmp1,*) 'repere : e r r o r  in kchoix : ', kchoix 
    340          CALL ctl_stop( ctmp1, ' kchoix must be equal to -1, 0 or 1 ' ) 
    341       ENDIF 
    342  
    343       ! 0.2 Initialization 
    344  
    345       zxx = gsinus(jpi/2,jpj/2)**2+gcosin(jpi/2,jpj/2)**2 
    346       IF( zxx <= 0.1 ) THEN 
    347          IF(lwp) WRITE(numout,*) 'repere : initialization ' 
    348          DO jj = 1, jpj 
    349             DO ji = 2, jpi 
    350                zlamdu(ji,jj) = glamu(ji,jj) - glamu(ji-1,jj) 
    351                zlamdu(ji,jj) = ASIN( SIN( rad*zlamdu(ji,jj) ) )/rad 
    352                zphiu (ji,jj) = gphiu(ji,jj) - gphiu(ji-1,jj) 
    353             END DO 
    354          END DO 
    355          DO jj = 2, jpj 
    356             zlamdv(:,jj) = glamv(:,jj)-glamv(:,jj-1) 
    357             zlamdv(:,jj) = ASIN(SIN(rad*zlamdv(:,jj)))/rad 
    358             zphiv (:,jj) = gphiv(:,jj)-gphiv(:,jj-1) 
    359          END DO 
    360           
    361          ! 0.3 Boudary conditions and periodicity 
    362           
    363          IF( nperio == 1 .OR.nperio == 4.OR.nperio == 6 ) THEN 
    364             DO jj = 1, jpj 
    365                zlamdu(1,jj) = zlamdu(jpim1,jj) 
    366                zphiu (1,jj) = zphiu (jpim1,jj) 
    367             END DO 
    368          ELSE 
    369             DO jj = 1, jpj 
    370                zlamdu(1,jj) = zlamdu(2,jj) 
    371                zphiu (1,jj) = zphiu (2,jj) 
    372             END DO 
    373          ENDIF 
    374           
    375          DO ji = 1, jpi 
    376             zlamdv(ji,1) = zlamdv(ji,2) 
    377             zphiv (ji,1) = zphiv (ji,2) 
    378          END DO 
    379           
    380          IF( nperio == 2 ) THEN 
    381             DO ji = 1, jpi 
    382                zphiv (ji,1) = zphiv (ji,3) 
    383             END DO 
    384          ENDIF 
    385           
    386          ! 0.4 gsinus gcosin 
    387           
    388 !CDIR NOVERRCHK 
    389          DO jj = 1, jpj 
    390 !CDIR NOVERRCHK 
    391             DO ji = 1, jpi 
    392                zcof1 = rad * ra * COS( rad * gphit(ji,jj) ) 
    393                zcof2 = rad * ra 
    394                zlamdu(ji,jj) = zlamdu(ji,jj) * zcof1 
    395                zlamdv(ji,jj) = zlamdv(ji,jj) * zcof1 
    396                zphiu (ji,jj) = zphiu (ji,jj) * zcof2 
    397                zphiv (ji,jj) = zphiv (ji,jj) * zcof2 
    398             END DO 
    399          END DO 
    400  
    401 !CDIR NOVERRCHK 
    402          DO jj = 1, jpj 
    403 !CDIR NOVERRCHK 
    404             DO ji = 1, jpi 
    405                ze1t = SQRT( zlamdu(ji,jj)*zlamdu(ji,jj) + zphiu(ji,jj)*zphiu(ji,jj) ) 
    406                ze2t = SQRT( zlamdv(ji,jj)*zlamdv(ji,jj) + zphiv(ji,jj)*zphiv(ji,jj) ) 
    407                gsinus(ji,jj) = 0.5*( zphiu(ji,jj)/ze1t - zlamdv(ji,jj)/ze2t ) 
    408                gcosin(ji,jj) = 0.5*( zphiv(ji,jj)/ze2t + zlamdu(ji,jj)/ze1t ) 
    409             END DO 
    410          END DO 
    411          CALL lbc_lnk( gsinus, 'T', -1. ) 
    412          CALL lbc_lnk( gcosin, 'T', -1. ) 
    413           
    414       ENDIF 
    415  
    416  
    417       ! 1. Change from geographic to tretched 
    418       ! ------------------------------------- 
    419        
    420       IF( kchoix == 1 ) THEN 
    421           px2(:,:) =  px1(:,:) * gcosin(:,:) + py1(:,:) * gsinus(:,:) 
    422           py2(:,:) = -px1(:,:) * gsinus(:,:) + py1(:,:) * gcosin(:,:) 
    423       ENDIF 
    424        
    425  
    426       ! 2. Change from tretched to geographic 
    427       ! ------------------------------------- 
    428        
    429       IF( kchoix == -1 ) THEN 
    430           px2(:,:) =  px1(:,:) * gcosin(:,:) - py1(:,:) * gsinus(:,:) 
    431           py2(:,:) =  px1(:,:) * gsinus(:,:) + py1(:,:) * gcosin(:,:) 
    432       ENDIF 
     430      CHARACTER(len=1), INTENT( IN ), OPTIONAL ::   cd_type      ! define the nature of pt2d array grid-points 
     431      ! 
     432      CHARACTER(len=1) ::   cl_type      ! define the nature of pt2d array grid-points (T point by default) 
     433      !!---------------------------------------------------------------------- 
     434 
     435      cl_type = 'T' 
     436      IF( PRESENT(cd_type) )   cl_type = cd_type 
     437         ! 
     438      SELECT CASE (kchoix) 
     439      CASE ( 1)      ! change from geographic to model grid. 
     440         px2(:,:) = rot_rep( px1, py1, cl_type, 'en->i' ) 
     441         py2(:,:) = rot_rep( px1, py1, cl_type, 'en->j' ) 
     442      CASE (-1)      ! change from model to geographic grid 
     443         px2(:,:) = rot_rep( px1, py1, cl_type, 'ij->e' ) 
     444         py2(:,:) = rot_rep( px1, py1, cl_type, 'ij->n' ) 
     445      CASE DEFAULT   ;   CALL ctl_stop( 'repere: Syntax Error in the definition of kchoix (1 OR -1' ) 
     446      END SELECT 
    433447       
    434448   END SUBROUTINE repere 
Note: See TracChangeset for help on using the changeset viewer.