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 5883 for branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90 – NEMO

Ignore:
Timestamp:
2015-11-13T08:01:08+01:00 (8 years ago)
Author:
gm
Message:

#1613: vvl by default: TRA/TRC remove optimization associated with linear free surface

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90

    r5836 r5883  
    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(:,:) 
     
    145114      CASE DEFAULT   ;   CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' ) 
    146115      END SELECT 
    147        
     116      ! 
    148117   END SUBROUTINE rot_rep 
    149118 
     
    155124      !! ** Purpose :   Compute angles between model grid lines and the North direction 
    156125      !! 
    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  
     126      !! ** Method  :   sinus and cosinus of the angle between the north-south axe  
     127      !!              and the j-direction at t, u, v and f-points 
     128      !!                dot and cross products are used to obtain cos and sin, resp. 
     129      !! 
     130      !! ** Action  : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf 
     131      !!---------------------------------------------------------------------- 
     132      INTEGER  ::   ji, jj   ! dummy loop indices 
     133      INTEGER  ::   ierr     ! local integer 
     134      REAL(wp) ::   zlam, zphi            ! local scalars 
     135      REAL(wp) ::   zlan, zphh            !   -      - 
     136      REAL(wp) ::   zxnpt, zynpt, znnpt   ! x,y components and norm of the vector: T point to North Pole 
     137      REAL(wp) ::   zxnpu, zynpu, znnpu   ! x,y components and norm of the vector: U point to North Pole 
     138      REAL(wp) ::   zxnpv, zynpv, znnpv   ! x,y components and norm of the vector: V point to North Pole 
     139      REAL(wp) ::   zxnpf, zynpf, znnpf   ! x,y components and norm of the vector: F point to North Pole 
     140      REAL(wp) ::   zxvvt, zyvvt, znvvt   ! x,y components and norm of the vector: between V points below and above a T point 
     141      REAL(wp) ::   zxffu, zyffu, znffu   ! x,y components and norm of the vector: between F points below and above a U point 
     142      REAL(wp) ::   zxffv, zyffv, znffv   ! x,y components and norm of the vector: between F points left  and right a V point 
     143      REAL(wp) ::   zxuuf, zyuuf, znuuf   ! x,y components and norm of the vector: between U points below and above a F point 
     144      !!---------------------------------------------------------------------- 
     145      ! 
    184146      ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj),   &  
    185147         &      gsinu(jpi,jpj), gcosu(jpi,jpj),   &  
     
    187149         &      gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr ) 
    188150      IF(lk_mpp)   CALL mpp_sum( ierr ) 
    189       IF( ierr /= 0 )   CALL ctl_stop('angle: unable to allocate arrays' ) 
    190  
     151      IF( ierr /= 0 )   CALL ctl_stop( 'angle: unable to allocate arrays' ) 
     152      ! 
    191153      ! ============================= ! 
    192154      ! Compute the cosinus and sinus ! 
    193155      ! ============================= ! 
    194156      ! (computation done on the north stereographic polar plane) 
    195  
     157      ! 
    196158      DO jj = 2, jpjm1 
    197159         DO ji = fs_2, jpi   ! vector opt. 
    198  
    199             ! north pole direction & modulous (at t-point) 
    200             zlam = glamt(ji,jj) 
     160            !                   
     161            zlam = glamt(ji,jj)     ! north pole direction & modulous (at t-point) 
    201162            zphi = gphit(ji,jj) 
    202163            zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    203164            zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    204165            znnpt = zxnpt*zxnpt + zynpt*zynpt 
    205  
    206             ! north pole direction & modulous (at u-point) 
    207             zlam = glamu(ji,jj) 
     166            ! 
     167            zlam = glamu(ji,jj)     ! north pole direction & modulous (at u-point) 
    208168            zphi = gphiu(ji,jj) 
    209169            zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    210170            zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    211171            znnpu = zxnpu*zxnpu + zynpu*zynpu 
    212  
    213             ! north pole direction & modulous (at v-point) 
    214             zlam = glamv(ji,jj) 
     172            ! 
     173            zlam = glamv(ji,jj)     ! north pole direction & modulous (at v-point) 
    215174            zphi = gphiv(ji,jj) 
    216175            zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    217176            zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    218177            znnpv = zxnpv*zxnpv + zynpv*zynpv 
    219  
    220             ! north pole direction & modulous (at f-point) 
    221             zlam = glamf(ji,jj) 
     178            ! 
     179            zlam = glamf(ji,jj)     ! north pole direction & modulous (at f-point) 
    222180            zphi = gphif(ji,jj) 
    223181            zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    224182            zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    225183            znnpf = zxnpf*zxnpf + zynpf*zynpf 
    226  
    227             ! j-direction: v-point segment direction (around t-point) 
    228             zlam = glamv(ji,jj  ) 
     184            ! 
     185            zlam = glamv(ji,jj  )   ! j-direction: v-point segment direction (around t-point) 
    229186            zphi = gphiv(ji,jj  ) 
    230187            zlan = glamv(ji,jj-1) 
     
    236193            znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  ) 
    237194            znvvt = MAX( znvvt, 1.e-14 ) 
    238  
    239             ! j-direction: f-point segment direction (around u-point) 
    240             zlam = glamf(ji,jj  ) 
     195            ! 
     196            zlam = glamf(ji,jj  )   ! j-direction: f-point segment direction (around u-point) 
    241197            zphi = gphif(ji,jj  ) 
    242198            zlan = glamf(ji,jj-1) 
     
    248204            znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu )  ) 
    249205            znffu = MAX( znffu, 1.e-14 ) 
    250  
    251             ! i-direction: f-point segment direction (around v-point) 
    252             zlam = glamf(ji  ,jj) 
     206            ! 
     207            zlam = glamf(ji  ,jj)   ! i-direction: f-point segment direction (around v-point) 
    253208            zphi = gphif(ji  ,jj) 
    254209            zlan = glamf(ji-1,jj) 
     
    260215            znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv )  ) 
    261216            znffv = MAX( znffv, 1.e-14 ) 
    262  
    263             ! j-direction: u-point segment direction (around f-point) 
    264             zlam = glamu(ji,jj+1) 
     217            ! 
     218            zlam = glamu(ji,jj+1)   ! j-direction: u-point segment direction (around f-point) 
    265219            zphi = gphiu(ji,jj+1) 
    266220            zlan = glamu(ji,jj  ) 
     
    272226            znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  ) 
    273227            znuuf = MAX( znuuf, 1.e-14 ) 
    274  
    275             ! cosinus and sinus using scalar and vectorial products 
     228            ! 
     229            !                       ! cosinus and sinus using dot and cross products 
    276230            gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt 
    277231            gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt 
    278  
     232            ! 
    279233            gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu 
    280234            gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu 
    281  
     235            ! 
    282236            gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf 
    283237            gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf 
    284  
    285             ! (caution, rotation of 90 degres) 
     238            ! 
    286239            gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv 
    287             gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv 
    288  
     240            gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv     ! (caution, rotation of 90 degres) 
     241            ! 
    289242         END DO 
    290243      END DO 
     
    318271      ! Lateral boundary conditions ! 
    319272      ! =========================== ! 
    320  
    321       ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 
     273      !           ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 
    322274      CALL lbc_lnk( gcost, 'T', -1. )   ;   CALL lbc_lnk( gsint, 'T', -1. ) 
    323275      CALL lbc_lnk( gcosu, 'U', -1. )   ;   CALL lbc_lnk( gsinu, 'U', -1. ) 
    324276      CALL lbc_lnk( gcosv, 'V', -1. )   ;   CALL lbc_lnk( gsinv, 'V', -1. ) 
    325277      CALL lbc_lnk( gcosf, 'F', -1. )   ;   CALL lbc_lnk( gsinf, 'F', -1. ) 
    326  
     278      ! 
    327279   END SUBROUTINE angle 
    328280 
    329281 
    330    SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid,     & 
    331                         pte, ptn ) 
     282   SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, pte, ptn ) 
    332283      !!---------------------------------------------------------------------- 
    333284      !!                    ***  ROUTINE geo2oce  *** 
     
    335286      !! ** Purpose : 
    336287      !! 
    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 
     288      !! ** Method  :   Change a vector from geocentric to east/north  
     289      !! 
    347290      !!---------------------------------------------------------------------- 
    348291      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::  pxx, pyy, pzz 
    349292      CHARACTER(len=1)            , INTENT(in   ) ::  cgrid 
    350293      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::  pte, ptn 
    351       !! 
     294      ! 
    352295      REAL(wp), PARAMETER :: rpi = 3.141592653e0 
    353296      REAL(wp), PARAMETER :: rad = rpi / 180.e0 
     
    355298      INTEGER ::   ierr   ! local integer 
    356299      !!---------------------------------------------------------------------- 
    357  
     300      ! 
    358301      IF( .NOT. ALLOCATED( gsinlon ) ) THEN 
    359302         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   & 
     
    361304         IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    362305         IF( ierr /= 0 )   CALL ctl_stop('geo2oce: unable to allocate arrays' ) 
     306      ENDIF 
     307      ! 
     308      SELECT CASE( cgrid) 
     309      CASE ( 'T' )    
     310         ig = 1 
     311         IF( .NOT. linit(ig) ) THEN  
     312            gsinlon(:,:,ig) = SIN( rad * glamt(:,:) ) 
     313            gcoslon(:,:,ig) = COS( rad * glamt(:,:) ) 
     314            gsinlat(:,:,ig) = SIN( rad * gphit(:,:) ) 
     315            gcoslat(:,:,ig) = COS( rad * gphit(:,:) ) 
     316            linit(ig) = .TRUE. 
     317         ENDIF 
     318      CASE ( 'U' )    
     319         ig = 2 
     320         IF( .NOT. linit(ig) ) THEN  
     321            gsinlon(:,:,ig) = SIN( rad * glamu(:,:) ) 
     322            gcoslon(:,:,ig) = COS( rad * glamu(:,:) ) 
     323            gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) ) 
     324            gcoslat(:,:,ig) = COS( rad * gphiu(:,:) ) 
     325            linit(ig) = .TRUE. 
     326         ENDIF 
     327      CASE ( 'V' )    
     328         ig = 3 
     329         IF( .NOT. linit(ig) ) THEN  
     330            gsinlon(:,:,ig) = SIN( rad * glamv(:,:) ) 
     331            gcoslon(:,:,ig) = COS( rad * glamv(:,:) ) 
     332            gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) ) 
     333            gcoslat(:,:,ig) = COS( rad * gphiv(:,:) ) 
     334            linit(ig) = .TRUE. 
     335         ENDIF 
     336      CASE ( 'F' )    
     337         ig = 4 
     338         IF( .NOT. linit(ig) ) THEN  
     339            gsinlon(:,:,ig) = SIN( rad * glamf(:,:) ) 
     340            gcoslon(:,:,ig) = COS( rad * glamf(:,:) ) 
     341            gsinlat(:,:,ig) = SIN( rad * gphif(:,:) ) 
     342            gcoslat(:,:,ig) = COS( rad * gphif(:,:) ) 
     343            linit(ig) = .TRUE. 
     344         ENDIF 
     345      CASE default    
     346         WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid 
     347         CALL ctl_stop( ctmp1 ) 
     348      END SELECT 
     349      ! 
     350      pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy 
     351      ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx    & 
     352         &  - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy    & 
     353         &  + gcoslat(:,:,ig) * pzz 
     354      ! 
     355   END SUBROUTINE geo2oce 
     356 
     357 
     358   SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz ) 
     359      !!---------------------------------------------------------------------- 
     360      !!                    ***  ROUTINE oce2geo  *** 
     361      !!       
     362      !! ** Purpose : 
     363      !! 
     364      !! ** Method  :   Change vector from east/north to geocentric 
     365      !! 
     366      !! History :     ! (A. Caubel)  oce2geo - Original code 
     367      !!---------------------------------------------------------------------- 
     368      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN    ) ::  pte, ptn 
     369      CHARACTER(len=1)            , INTENT( IN    ) ::  cgrid 
     370      REAL(wp), DIMENSION(jpi,jpj), INTENT(   OUT ) ::  pxx , pyy , pzz 
     371      !! 
     372      REAL(wp), PARAMETER :: rpi = 3.141592653E0 
     373      REAL(wp), PARAMETER :: rad = rpi / 180.e0 
     374      INTEGER ::   ig     ! 
     375      INTEGER ::   ierr   ! local integer 
     376      !!---------------------------------------------------------------------- 
     377 
     378      IF( .NOT. ALLOCATED( gsinlon ) ) THEN 
     379         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   & 
     380            &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 
     381         IF( lk_mpp    )   CALL mpp_sum( ierr ) 
     382         IF( ierr /= 0 )   CALL ctl_stop('oce2geo: unable to allocate arrays' ) 
    363383      ENDIF 
    364384 
     
    404424            CALL ctl_stop( ctmp1 ) 
    405425      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        
     426      ! 
     427      pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn  
     428      pyy =   gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn 
     429      pzz =   gcoslat(:,:,ig) * ptn 
     430      ! 
    493431   END SUBROUTINE oce2geo 
    494432 
    495433 
    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 ) 
     434   SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv ) 
    541435      !!---------------------------------------------------------------------- 
    542436      !!                  ***  ROUTINE obs_rot  *** 
     
    546440      !!                current at observation points 
    547441      !! 
    548       !! History : 
    549       !!   9.2  !  09-02  (K. Mogensen) 
     442      !! History :  9.2  !  09-02  (K. Mogensen) 
    550443      !!---------------------------------------------------------------------- 
    551444      REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT )::   psinu, pcosu, psinv, pcosv   ! copy of data 
    552445      !!---------------------------------------------------------------------- 
    553  
     446      ! 
    554447      ! Initialization of gsin* and gcos* at first call 
    555448      ! ----------------------------------------------- 
    556  
    557449      IF( lmust_init ) THEN 
    558450         IF(lwp) WRITE(numout,*) 
    559451         IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched' 
    560452         IF(lwp) WRITE(numout,*) ' ~~~~~~~   coordinate transformation' 
    561  
    562453         CALL angle       ! initialization of the transformation 
    563454         lmust_init = .FALSE. 
    564  
    565455      ENDIF 
    566  
     456      ! 
    567457      psinu(:,:) = gsinu(:,:) 
    568458      pcosu(:,:) = gcosu(:,:) 
    569459      psinv(:,:) = gsinv(:,:) 
    570460      pcosv(:,:) = gcosv(:,:) 
    571  
     461      ! 
    572462   END SUBROUTINE obs_rot 
    573463 
Note: See TracChangeset for help on using the changeset viewer.