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 7351 for branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 – NEMO

Ignore:
Timestamp:
2016-11-28T17:04:10+01:00 (7 years ago)
Author:
emanuelaclementi
Message:

ticket #1805 step 3: /2016/dev_INGV_UKMO_2016 aligned to the trunk at revision 7161

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90

    r5836 r7351  
    55   !!====================================================================== 
    66   !! History :  OPA  !  07-1996  (O. Marti)  Original code 
    7    !!   NEMO     1.0  !  02-2008  (G. Madec)  F90: Free form 
    8    !!            3.0  !   
     7   !!   NEMO     1.0  !  06-2006  (G. Madec )  Free form, F90 + opt. 
     8   !!                 !  04-2007  (S. Masson)  angle: Add T, F points and bugfix in cos lateral boundary 
     9   !!            3.0  !  07-2008  (G. Madec)  geo2oce suppress lon/lat agruments 
     10   !!            3.7  !  11-2015  (G. Madec)  remove the unused repere and repcmo routines 
    911   !!---------------------------------------------------------------------- 
    1012 
    1113   !!---------------------------------------------------------------------- 
    12    !!   repcmo      :  
    13    !!   angle       : 
    14    !!   geo2oce     : 
    15    !!   repere      :   old routine suppress it ??? 
     14   !!   rot_rep       : Rotate the Repere: geographic grid <==> stretched coordinates grid 
     15   !!   angle         : 
     16   !!   geo2oce       : 
     17   !!   oce2geo       : 
    1618   !!---------------------------------------------------------------------- 
    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 
     19   USE dom_oce        ! mesh and scale factors 
     20   USE phycst         ! physical constants 
     21   ! 
     22   USE in_out_manager ! I/O manager 
     23   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     24   USE lib_mpp        ! MPP library 
    2225 
    2326   IMPLICIT NONE 
    2427   PRIVATE 
    2528 
    26    PUBLIC   rot_rep, repcmo, repere, geo2oce, oce2geo   ! only rot_rep should be used 
    27                                              ! repcmo and repere are keep only for compatibility. 
    28                                              ! they are only a useless overlay of rot_rep 
    29  
    30    PUBLIC   obs_rot 
    31  
    32    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   & 
    33       gsint, gcost,   &  ! cos/sin between model grid lines and NP direction at T point 
    34       gsinu, gcosu,   &  ! cos/sin between model grid lines and NP direction at U point 
    35       gsinv, gcosv,   &  ! cos/sin between model grid lines and NP direction at V point 
    36       gsinf, gcosf       ! cos/sin between model grid lines and NP direction at F point 
     29   PUBLIC   rot_rep   ! called in sbccpl, fldread, and cyclone 
     30   PUBLIC   geo2oce   ! called in sbccpl 
     31   PUBLIC   oce2geo   ! called in sbccpl 
     32   PUBLIC   obs_rot   ! called in obs_rot_vel and obs_write 
     33 
     34   !                                         ! cos/sin between model grid lines and NP direction 
     35   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsint, gcost   ! at T point 
     36   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinu, gcosu   ! at U point 
     37   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinv, gcosv   ! at V point 
     38   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gsinf, gcosf   ! at F point 
    3739 
    3840   LOGICAL ,              SAVE, DIMENSION(4)     ::   linit = .FALSE. 
    3941   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gsinlon, gcoslon, gsinlat, gcoslat 
    4042 
    41    LOGICAL ::   lmust_init = .TRUE.        !: used to initialize the cos/sin variables (se above) 
     43   LOGICAL ::   lmust_init = .TRUE.        !: used to initialize the cos/sin variables (see above) 
    4244 
    4345   !! * Substitutions 
     
    5052CONTAINS 
    5153 
    52    SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1,   & 
    53                        px2 , py2 ) 
    54       !!---------------------------------------------------------------------- 
    55       !!                  ***  ROUTINE repcmo  *** 
    56       !! 
    57       !! ** Purpose :   Change vector componantes from a geographic grid to a 
    58       !!      stretched coordinates grid. 
    59       !! 
    60       !! ** Method  :   Initialization of arrays at the first call. 
    61       !! 
    62       !! ** Action  : - px2 : first  componante (defined at u point) 
    63       !!              - py2 : second componante (defined at v point) 
    64       !!---------------------------------------------------------------------- 
    65       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pxu1, pyu1   ! geographic vector componantes at u-point 
    66       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pxv1, pyv1   ! geographic vector componantes at v-point 
    67       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   px2          ! i-componante (defined at u-point) 
    68       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   py2          ! j-componante (defined at v-point) 
    69       !!---------------------------------------------------------------------- 
    70        
    71       ! Change from geographic to stretched coordinate 
    72       ! ---------------------------------------------- 
    73       CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 ) 
    74       CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 ) 
    75        
    76    END SUBROUTINE repcmo 
    77  
    78  
    7954   SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot ) 
    8055      !!---------------------------------------------------------------------- 
     
    8358      !! ** Purpose :   Rotate the Repere: Change vector componantes between 
    8459      !!                geographic grid <--> stretched coordinates grid. 
    85       !! 
    86       !! History : 
    87       !!   9.2  !  07-04  (S. Masson)   
    88       !!                  (O. Marti ) Original code (repere and repcmo) 
    89       !!---------------------------------------------------------------------- 
    90       REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) ::   pxin, pyin   ! vector componantes 
    91       CHARACTER(len=1),             INTENT( IN ) ::   cd_type      ! define the nature of pt2d array grid-points 
    92       CHARACTER(len=5),             INTENT( IN ) ::   cdtodo       ! specify the work to do: 
    93       !!                                                           ! 'en->i' east-north componantes to model i componante 
    94       !!                                                           ! 'en->j' east-north componantes to model j componante 
    95       !!                                                           ! 'ij->e' model i-j componantes to east componante 
    96       !!                                                           ! 'ij->n' model i-j componantes to east componante 
    97       REAL(wp), DIMENSION(jpi,jpj), INTENT(out) ::   prot       
    98       !!---------------------------------------------------------------------- 
    99  
    100       ! Initialization of gsin* and gcos* at first call 
    101       ! ----------------------------------------------- 
    102  
    103       IF( lmust_init ) THEN 
     60      !!---------------------------------------------------------------------- 
     61      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pxin, pyin   ! vector componantes 
     62      CHARACTER(len=1),             INTENT(in   ) ::   cd_type      ! define the nature of pt2d array grid-points 
     63      CHARACTER(len=5),             INTENT(in   ) ::   cdtodo       ! type of transpormation: 
     64      !                                                             ! 'en->i' = east-north to i-component 
     65      !                                                             ! 'en->j' = east-north to j-component 
     66      !                                                             ! 'ij->e' = (i,j) components to east 
     67      !                                                             ! 'ij->n' = (i,j) components to north 
     68      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   prot       
     69      !!---------------------------------------------------------------------- 
     70      ! 
     71      IF( lmust_init ) THEN      ! at 1st call only: set  gsin. & gcos. 
    10472         IF(lwp) WRITE(numout,*) 
    105          IF(lwp) WRITE(numout,*) ' rot_rep : geographic <--> stretched' 
    106          IF(lwp) WRITE(numout,*) ' ~~~~~    coordinate transformation' 
     73         IF(lwp) WRITE(numout,*) ' rot_rep: coordinate transformation : geographic <==> model (i,j)-components' 
     74         IF(lwp) WRITE(numout,*) ' ~~~~~~~~    ' 
    10775         ! 
    10876         CALL angle       ! initialization of the transformation 
    10977         lmust_init = .FALSE. 
    11078      ENDIF 
    111        
    112       SELECT CASE (cdtodo) 
    113       CASE ('en->i')      ! 'en->i' est-north componantes to model i componante 
     79      ! 
     80      SELECT CASE( cdtodo )      ! type of rotation 
     81      ! 
     82      CASE( 'en->i' )                  ! east-north to i-component 
    11483         SELECT CASE (cd_type) 
    11584         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:) 
     
    11988         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 
    12089         END SELECT 
    121       CASE ('en->j')      ! 'en->j' est-north componantes to model j componante 
     90      CASE ('en->j')                   ! east-north to j-component 
    12291         SELECT CASE (cd_type) 
    12392         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:) 
     
    12796         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 
    12897         END SELECT 
    129       CASE ('ij->e')      ! 'ij->e' model i-j componantes to est componante 
     98      CASE ('ij->e')                   ! (i,j)-components to east 
    13099         SELECT CASE (cd_type) 
    131100         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:) 
     
    135104         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 
    136105         END SELECT 
    137       CASE ('ij->n')      ! 'ij->n' model i-j componantes to est componante 
     106      CASE ('ij->n')                   ! (i,j)-components to north  
    138107         SELECT CASE (cd_type) 
    139108         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:) 
     
    144113         END SELECT 
    145114      CASE DEFAULT   ;   CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' ) 
     115      ! 
    146116      END SELECT 
    147        
     117      ! 
    148118   END SUBROUTINE rot_rep 
    149119 
     
    155125      !! ** Purpose :   Compute angles between model grid lines and the North direction 
    156126      !! 
    157       !! ** Method  : 
    158       !! 
    159       !! ** Action  :   Compute (gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf) arrays: 
    160       !!      sinus and cosinus of the angle between the north-south axe and the  
    161       !!      j-direction at t, u, v and f-points 
    162       !! 
    163       !! History : 
    164       !!   7.0  !  96-07  (O. Marti )  Original code 
    165       !!   8.0  !  98-06  (G. Madec ) 
    166       !!   8.5  !  98-06  (G. Madec )  Free form, F90 + opt. 
    167       !!   9.2  !  07-04  (S. Masson)  Add T, F points and bugfix in cos lateral boundary 
    168       !!---------------------------------------------------------------------- 
    169       INTEGER ::   ji, jj   ! dummy loop indices 
    170       INTEGER ::   ierr     ! local integer 
    171       REAL(wp) ::   & 
    172          zlam, zphi,            &  ! temporary scalars 
    173          zlan, zphh,            &  !    "         " 
    174          zxnpt, zynpt, znnpt,   &  ! x,y components and norm of the vector: T point to North Pole 
    175          zxnpu, zynpu, znnpu,   &  ! x,y components and norm of the vector: U point to North Pole 
    176          zxnpv, zynpv, znnpv,   &  ! x,y components and norm of the vector: V point to North Pole 
    177          zxnpf, zynpf, znnpf,   &  ! x,y components and norm of the vector: F point to North Pole 
    178          zxvvt, zyvvt, znvvt,   &  ! x,y components and norm of the vector: between V points below and above a T point 
    179          zxffu, zyffu, znffu,   &  ! x,y components and norm of the vector: between F points below and above a U point 
    180          zxffv, zyffv, znffv,   &  ! x,y components and norm of the vector: between F points left  and right a V point 
    181          zxuuf, zyuuf, znuuf       ! x,y components and norm of the vector: between U points below and above a F point 
    182       !!---------------------------------------------------------------------- 
    183  
     127      !! ** Method  :   sinus and cosinus of the angle between the north-south axe  
     128      !!              and the j-direction at t, u, v and f-points 
     129      !!                dot and cross products are used to obtain cos and sin, resp. 
     130      !! 
     131      !! ** Action  : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf 
     132      !!---------------------------------------------------------------------- 
     133      INTEGER  ::   ji, jj   ! dummy loop indices 
     134      INTEGER  ::   ierr     ! local integer 
     135      REAL(wp) ::   zlam, zphi            ! local scalars 
     136      REAL(wp) ::   zlan, zphh            !   -      - 
     137      REAL(wp) ::   zxnpt, zynpt, znnpt   ! x,y components and norm of the vector: T point to North Pole 
     138      REAL(wp) ::   zxnpu, zynpu, znnpu   ! x,y components and norm of the vector: U point to North Pole 
     139      REAL(wp) ::   zxnpv, zynpv, znnpv   ! x,y components and norm of the vector: V point to North Pole 
     140      REAL(wp) ::   zxnpf, zynpf, znnpf   ! x,y components and norm of the vector: F point to North Pole 
     141      REAL(wp) ::   zxvvt, zyvvt, znvvt   ! x,y components and norm of the vector: between V points below and above a T point 
     142      REAL(wp) ::   zxffu, zyffu, znffu   ! x,y components and norm of the vector: between F points below and above a U point 
     143      REAL(wp) ::   zxffv, zyffv, znffv   ! x,y components and norm of the vector: between F points left  and right a V point 
     144      REAL(wp) ::   zxuuf, zyuuf, znuuf   ! x,y components and norm of the vector: between U points below and above a F point 
     145      !!---------------------------------------------------------------------- 
     146      ! 
    184147      ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj),   &  
    185148         &      gsinu(jpi,jpj), gcosu(jpi,jpj),   &  
     
    187150         &      gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr ) 
    188151      IF(lk_mpp)   CALL mpp_sum( ierr ) 
    189       IF( ierr /= 0 )   CALL ctl_stop('angle: unable to allocate arrays' ) 
    190  
     152      IF( ierr /= 0 )   CALL ctl_stop( 'angle: unable to allocate arrays' ) 
     153      ! 
    191154      ! ============================= ! 
    192155      ! Compute the cosinus and sinus ! 
    193156      ! ============================= ! 
    194157      ! (computation done on the north stereographic polar plane) 
    195  
     158      ! 
    196159      DO jj = 2, jpjm1 
    197160         DO ji = fs_2, jpi   ! vector opt. 
    198  
    199             ! north pole direction & modulous (at t-point) 
    200             zlam = glamt(ji,jj) 
     161            !                   
     162            zlam = glamt(ji,jj)     ! north pole direction & modulous (at t-point) 
    201163            zphi = gphit(ji,jj) 
    202164            zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    203165            zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    204166            znnpt = zxnpt*zxnpt + zynpt*zynpt 
    205  
    206             ! north pole direction & modulous (at u-point) 
    207             zlam = glamu(ji,jj) 
     167            ! 
     168            zlam = glamu(ji,jj)     ! north pole direction & modulous (at u-point) 
    208169            zphi = gphiu(ji,jj) 
    209170            zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    210171            zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    211172            znnpu = zxnpu*zxnpu + zynpu*zynpu 
    212  
    213             ! north pole direction & modulous (at v-point) 
    214             zlam = glamv(ji,jj) 
     173            ! 
     174            zlam = glamv(ji,jj)     ! north pole direction & modulous (at v-point) 
    215175            zphi = gphiv(ji,jj) 
    216176            zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    217177            zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    218178            znnpv = zxnpv*zxnpv + zynpv*zynpv 
    219  
    220             ! north pole direction & modulous (at f-point) 
    221             zlam = glamf(ji,jj) 
     179            ! 
     180            zlam = glamf(ji,jj)     ! north pole direction & modulous (at f-point) 
    222181            zphi = gphif(ji,jj) 
    223182            zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    224183            zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    225184            znnpf = zxnpf*zxnpf + zynpf*zynpf 
    226  
    227             ! j-direction: v-point segment direction (around t-point) 
    228             zlam = glamv(ji,jj  ) 
     185            ! 
     186            zlam = glamv(ji,jj  )   ! j-direction: v-point segment direction (around t-point) 
    229187            zphi = gphiv(ji,jj  ) 
    230188            zlan = glamv(ji,jj-1) 
     
    236194            znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  ) 
    237195            znvvt = MAX( znvvt, 1.e-14 ) 
    238  
    239             ! j-direction: f-point segment direction (around u-point) 
    240             zlam = glamf(ji,jj  ) 
     196            ! 
     197            zlam = glamf(ji,jj  )   ! j-direction: f-point segment direction (around u-point) 
    241198            zphi = gphif(ji,jj  ) 
    242199            zlan = glamf(ji,jj-1) 
     
    248205            znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu )  ) 
    249206            znffu = MAX( znffu, 1.e-14 ) 
    250  
    251             ! i-direction: f-point segment direction (around v-point) 
    252             zlam = glamf(ji  ,jj) 
     207            ! 
     208            zlam = glamf(ji  ,jj)   ! i-direction: f-point segment direction (around v-point) 
    253209            zphi = gphif(ji  ,jj) 
    254210            zlan = glamf(ji-1,jj) 
     
    260216            znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv )  ) 
    261217            znffv = MAX( znffv, 1.e-14 ) 
    262  
    263             ! j-direction: u-point segment direction (around f-point) 
    264             zlam = glamu(ji,jj+1) 
     218            ! 
     219            zlam = glamu(ji,jj+1)   ! j-direction: u-point segment direction (around f-point) 
    265220            zphi = gphiu(ji,jj+1) 
    266221            zlan = glamu(ji,jj  ) 
     
    272227            znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  ) 
    273228            znuuf = MAX( znuuf, 1.e-14 ) 
    274  
    275             ! cosinus and sinus using scalar and vectorial products 
     229            ! 
     230            !                       ! cosinus and sinus using dot and cross products 
    276231            gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt 
    277232            gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt 
    278  
     233            ! 
    279234            gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu 
    280235            gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu 
    281  
     236            ! 
    282237            gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf 
    283238            gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf 
    284  
    285             ! (caution, rotation of 90 degres) 
     239            ! 
    286240            gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv 
    287             gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv 
    288  
     241            gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv     ! (caution, rotation of 90 degres) 
     242            ! 
    289243         END DO 
    290244      END DO 
     
    318272      ! Lateral boundary conditions ! 
    319273      ! =========================== ! 
    320  
    321       ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 
     274      !           ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 
    322275      CALL lbc_lnk( gcost, 'T', -1. )   ;   CALL lbc_lnk( gsint, 'T', -1. ) 
    323276      CALL lbc_lnk( gcosu, 'U', -1. )   ;   CALL lbc_lnk( gsinu, 'U', -1. ) 
    324277      CALL lbc_lnk( gcosv, 'V', -1. )   ;   CALL lbc_lnk( gsinv, 'V', -1. ) 
    325278      CALL lbc_lnk( gcosf, 'F', -1. )   ;   CALL lbc_lnk( gsinf, 'F', -1. ) 
    326  
     279      ! 
    327280   END SUBROUTINE angle 
    328281 
    329282 
    330    SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid,     & 
    331                         pte, ptn ) 
     283   SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, pte, ptn ) 
    332284      !!---------------------------------------------------------------------- 
    333285      !!                    ***  ROUTINE geo2oce  *** 
     
    335287      !! ** Purpose : 
    336288      !! 
    337       !! ** Method  :   Change wind stress from geocentric to east/north 
    338       !! 
    339       !! History : 
    340       !!        !         (O. Marti)  Original code 
    341       !!        !  91-03  (G. Madec) 
    342       !!        !  92-07  (M. Imbard) 
    343       !!        !  99-11  (M. Imbard) NetCDF format with IOIPSL 
    344       !!        !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb 
    345       !!   8.5  !  02-06  (G. Madec)  F90: Free form 
    346       !!   3.0  !  07-08  (G. Madec)  geo2oce suppress lon/lat agruments 
     289      !! ** Method  :   Change a vector from geocentric to east/north  
     290      !! 
    347291      !!---------------------------------------------------------------------- 
    348292      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::  pxx, pyy, pzz 
    349293      CHARACTER(len=1)            , INTENT(in   ) ::  cgrid 
    350294      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::  pte, ptn 
    351       !! 
     295      ! 
    352296      REAL(wp), PARAMETER :: rpi = 3.141592653e0 
    353297      REAL(wp), PARAMETER :: rad = rpi / 180.e0 
     
    355299      INTEGER ::   ierr   ! local integer 
    356300      !!---------------------------------------------------------------------- 
    357  
     301      ! 
    358302      IF( .NOT. ALLOCATED( gsinlon ) ) THEN 
    359303         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   & 
     
    361305         IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    362306         IF( ierr /= 0 )   CALL ctl_stop('geo2oce: unable to allocate arrays' ) 
     307      ENDIF 
     308      ! 
     309      SELECT CASE( cgrid) 
     310      CASE ( 'T' )    
     311         ig = 1 
     312         IF( .NOT. linit(ig) ) THEN  
     313            gsinlon(:,:,ig) = SIN( rad * glamt(:,:) ) 
     314            gcoslon(:,:,ig) = COS( rad * glamt(:,:) ) 
     315            gsinlat(:,:,ig) = SIN( rad * gphit(:,:) ) 
     316            gcoslat(:,:,ig) = COS( rad * gphit(:,:) ) 
     317            linit(ig) = .TRUE. 
     318         ENDIF 
     319      CASE ( 'U' )    
     320         ig = 2 
     321         IF( .NOT. linit(ig) ) THEN  
     322            gsinlon(:,:,ig) = SIN( rad * glamu(:,:) ) 
     323            gcoslon(:,:,ig) = COS( rad * glamu(:,:) ) 
     324            gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) ) 
     325            gcoslat(:,:,ig) = COS( rad * gphiu(:,:) ) 
     326            linit(ig) = .TRUE. 
     327         ENDIF 
     328      CASE ( 'V' )    
     329         ig = 3 
     330         IF( .NOT. linit(ig) ) THEN  
     331            gsinlon(:,:,ig) = SIN( rad * glamv(:,:) ) 
     332            gcoslon(:,:,ig) = COS( rad * glamv(:,:) ) 
     333            gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) ) 
     334            gcoslat(:,:,ig) = COS( rad * gphiv(:,:) ) 
     335            linit(ig) = .TRUE. 
     336         ENDIF 
     337      CASE ( 'F' )    
     338         ig = 4 
     339         IF( .NOT. linit(ig) ) THEN  
     340            gsinlon(:,:,ig) = SIN( rad * glamf(:,:) ) 
     341            gcoslon(:,:,ig) = COS( rad * glamf(:,:) ) 
     342            gsinlat(:,:,ig) = SIN( rad * gphif(:,:) ) 
     343            gcoslat(:,:,ig) = COS( rad * gphif(:,:) ) 
     344            linit(ig) = .TRUE. 
     345         ENDIF 
     346      CASE default    
     347         WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid 
     348         CALL ctl_stop( ctmp1 ) 
     349      END SELECT 
     350      ! 
     351      pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy 
     352      ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx    & 
     353         &  - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy    & 
     354         &  + gcoslat(:,:,ig) * pzz 
     355      ! 
     356   END SUBROUTINE geo2oce 
     357 
     358 
     359   SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz ) 
     360      !!---------------------------------------------------------------------- 
     361      !!                    ***  ROUTINE oce2geo  *** 
     362      !!       
     363      !! ** Purpose : 
     364      !! 
     365      !! ** Method  :   Change vector from east/north to geocentric 
     366      !! 
     367      !! History :     ! (A. Caubel)  oce2geo - Original code 
     368      !!---------------------------------------------------------------------- 
     369      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN    ) ::  pte, ptn 
     370      CHARACTER(len=1)            , INTENT( IN    ) ::  cgrid 
     371      REAL(wp), DIMENSION(jpi,jpj), INTENT(   OUT ) ::  pxx , pyy , pzz 
     372      !! 
     373      REAL(wp), PARAMETER :: rpi = 3.141592653E0 
     374      REAL(wp), PARAMETER :: rad = rpi / 180.e0 
     375      INTEGER ::   ig     ! 
     376      INTEGER ::   ierr   ! local integer 
     377      !!---------------------------------------------------------------------- 
     378 
     379      IF( .NOT. ALLOCATED( gsinlon ) ) THEN 
     380         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   & 
     381            &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 
     382         IF( lk_mpp    )   CALL mpp_sum( ierr ) 
     383         IF( ierr /= 0 )   CALL ctl_stop('oce2geo: unable to allocate arrays' ) 
    363384      ENDIF 
    364385 
     
    404425            CALL ctl_stop( ctmp1 ) 
    405426      END SELECT 
    406        
    407       pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy 
    408       ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx    & 
    409             - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy    & 
    410             + gcoslat(:,:,ig) * pzz 
    411 !!$   ptv =   gcoslon(:,:,ig) * gcoslat(:,:,ig) * pxx    & 
    412 !!$         + gsinlon(:,:,ig) * gcoslat(:,:,ig) * pyy    & 
    413 !!$         + gsinlat(:,:,ig) * pzz 
    414       ! 
    415    END SUBROUTINE geo2oce 
    416  
    417    SUBROUTINE oce2geo ( pte, ptn, cgrid,     & 
    418                         pxx , pyy , pzz ) 
    419       !!---------------------------------------------------------------------- 
    420       !!                    ***  ROUTINE oce2geo  *** 
    421       !!       
    422       !! ** Purpose : 
    423       !! 
    424       !! ** Method  :   Change vector from east/north to geocentric 
    425       !! 
    426       !! History : 
    427       !!        !         (A. Caubel)  oce2geo - Original code 
    428       !!---------------------------------------------------------------------- 
    429       REAL(wp), DIMENSION(jpi,jpj), INTENT( IN    ) ::  pte, ptn 
    430       CHARACTER(len=1)            , INTENT( IN    ) ::  cgrid 
    431       REAL(wp), DIMENSION(jpi,jpj), INTENT(   OUT ) ::  pxx , pyy , pzz 
    432       !! 
    433       REAL(wp), PARAMETER :: rpi = 3.141592653E0 
    434       REAL(wp), PARAMETER :: rad = rpi / 180.e0 
    435       INTEGER ::   ig     ! 
    436       INTEGER ::   ierr   ! local integer 
    437       !!---------------------------------------------------------------------- 
    438  
    439       IF( .NOT. ALLOCATED( gsinlon ) ) THEN 
    440          ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   & 
    441             &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 
    442          IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    443          IF( ierr /= 0 )   CALL ctl_stop('oce2geo: unable to allocate arrays' ) 
    444       ENDIF 
    445  
    446       SELECT CASE( cgrid) 
    447          CASE ( 'T' )    
    448             ig = 1 
    449             IF( .NOT. linit(ig) ) THEN  
    450                gsinlon(:,:,ig) = SIN( rad * glamt(:,:) ) 
    451                gcoslon(:,:,ig) = COS( rad * glamt(:,:) ) 
    452                gsinlat(:,:,ig) = SIN( rad * gphit(:,:) ) 
    453                gcoslat(:,:,ig) = COS( rad * gphit(:,:) ) 
    454                linit(ig) = .TRUE. 
    455             ENDIF 
    456          CASE ( 'U' )    
    457             ig = 2 
    458             IF( .NOT. linit(ig) ) THEN  
    459                gsinlon(:,:,ig) = SIN( rad * glamu(:,:) ) 
    460                gcoslon(:,:,ig) = COS( rad * glamu(:,:) ) 
    461                gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) ) 
    462                gcoslat(:,:,ig) = COS( rad * gphiu(:,:) ) 
    463                linit(ig) = .TRUE. 
    464             ENDIF 
    465          CASE ( 'V' )    
    466             ig = 3 
    467             IF( .NOT. linit(ig) ) THEN  
    468                gsinlon(:,:,ig) = SIN( rad * glamv(:,:) ) 
    469                gcoslon(:,:,ig) = COS( rad * glamv(:,:) ) 
    470                gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) ) 
    471                gcoslat(:,:,ig) = COS( rad * gphiv(:,:) ) 
    472                linit(ig) = .TRUE. 
    473             ENDIF 
    474          CASE ( 'F' )    
    475             ig = 4 
    476             IF( .NOT. linit(ig) ) THEN  
    477                gsinlon(:,:,ig) = SIN( rad * glamf(:,:) ) 
    478                gcoslon(:,:,ig) = COS( rad * glamf(:,:) ) 
    479                gsinlat(:,:,ig) = SIN( rad * gphif(:,:) ) 
    480                gcoslat(:,:,ig) = COS( rad * gphif(:,:) ) 
    481                linit(ig) = .TRUE. 
    482             ENDIF 
    483          CASE default    
    484             WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid 
    485             CALL ctl_stop( ctmp1 ) 
    486       END SELECT 
    487  
    488        pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn  
    489        pyy =   gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn 
    490        pzz =   gcoslat(:,:,ig) * ptn 
    491  
    492        
     427      ! 
     428      pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn  
     429      pyy =   gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn 
     430      pzz =   gcoslat(:,:,ig) * ptn 
     431      ! 
    493432   END SUBROUTINE oce2geo 
    494433 
    495434 
    496    SUBROUTINE repere ( px1, py1, px2, py2, kchoix, cd_type ) 
    497       !!---------------------------------------------------------------------- 
    498       !!                 ***  ROUTINE repere  *** 
    499       !!         
    500       !! ** Purpose :   Change vector componantes between a geopgraphic grid  
    501       !!      and a stretched coordinates grid. 
    502       !! 
    503       !! ** Method  :    
    504       !! 
    505       !! ** Action  : 
    506       !! 
    507       !! History : 
    508       !!        !  89-03  (O. Marti)  original code 
    509       !!        !  92-02  (M. Imbard) 
    510       !!        !  93-03  (M. Guyon)  symetrical conditions 
    511       !!        !  98-05  (B. Blanke) 
    512       !!   8.5  !  02-08  (G. Madec)  F90: Free form 
    513       !!---------------------------------------------------------------------- 
    514       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   px1, py1   ! two horizontal components to be rotated 
    515       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   px2, py2   ! the two horizontal components in the model repere 
    516       INTEGER , INTENT(in   )                     ::   kchoix     ! type of transformation 
    517       !                                                           ! = 1 change from geographic to model grid. 
    518       !                                                           ! =-1 change from model to geographic grid 
    519       CHARACTER(len=1), INTENT(in   ), OPTIONAL   ::   cd_type    ! define the nature of pt2d array grid-points 
    520       ! 
    521       CHARACTER(len=1) ::   cl_type      ! define the nature of pt2d array grid-points (T point by default) 
    522       !!---------------------------------------------------------------------- 
    523  
    524       cl_type = 'T' 
    525       IF( PRESENT(cd_type) )   cl_type = cd_type 
    526          ! 
    527       SELECT CASE (kchoix) 
    528       CASE ( 1)      ! change from geographic to model grid. 
    529          CALL rot_rep( px1, py1, cl_type, 'en->i', px2 ) 
    530          CALL rot_rep( px1, py1, cl_type, 'en->j', py2 ) 
    531       CASE (-1)      ! change from model to geographic grid 
    532          CALL rot_rep( px1, py1, cl_type, 'ij->e', px2 ) 
    533          CALL rot_rep( px1, py1, cl_type, 'ij->n', py2 ) 
    534       CASE DEFAULT   ;   CALL ctl_stop( 'repere: Syntax Error in the definition of kchoix (1 OR -1' ) 
    535       END SELECT 
    536        
    537    END SUBROUTINE repere 
    538  
    539  
    540    SUBROUTINE obs_rot ( psinu, pcosu, psinv, pcosv ) 
     435   SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv ) 
    541436      !!---------------------------------------------------------------------- 
    542437      !!                  ***  ROUTINE obs_rot  *** 
     
    546441      !!                current at observation points 
    547442      !! 
    548       !! History : 
    549       !!   9.2  !  09-02  (K. Mogensen) 
     443      !! History :  9.2  !  09-02  (K. Mogensen) 
    550444      !!---------------------------------------------------------------------- 
    551445      REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT )::   psinu, pcosu, psinv, pcosv   ! copy of data 
    552446      !!---------------------------------------------------------------------- 
    553  
     447      ! 
    554448      ! Initialization of gsin* and gcos* at first call 
    555449      ! ----------------------------------------------- 
    556  
    557450      IF( lmust_init ) THEN 
    558451         IF(lwp) WRITE(numout,*) 
    559452         IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched' 
    560453         IF(lwp) WRITE(numout,*) ' ~~~~~~~   coordinate transformation' 
    561  
    562454         CALL angle       ! initialization of the transformation 
    563455         lmust_init = .FALSE. 
    564  
    565456      ENDIF 
    566  
     457      ! 
    567458      psinu(:,:) = gsinu(:,:) 
    568459      pcosu(:,:) = gcosu(:,:) 
    569460      psinv(:,:) = gsinv(:,:) 
    570461      pcosv(:,:) = gcosv(:,:) 
    571  
     462      ! 
    572463   END SUBROUTINE obs_rot 
    573464 
Note: See TracChangeset for help on using the changeset viewer.