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 6041 for branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 – NEMO

Ignore:
Timestamp:
2015-12-14T10:06:06+01:00 (8 years ago)
Author:
timgraham
Message:

Merged head of trunk into branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r5656 r6041  
    1414   !!                            use of parameters in par_CONFIG-Rxx.h90, not in namelist 
    1515   !!             -   ! 2004-05  (A. Koch-Larrouy) Add Gyre configuration  
    16    !!            4.0  ! 2011-02  (G. Madec) add cell surface (e1e2t) 
     16   !!            3.7  ! 2015-09  (G. Madec, S. Flavoni) add cell surface and their inverse 
     17   !!                                       add optional read of e1e2u & e1e2v 
    1718   !!---------------------------------------------------------------------- 
    1819 
     
    2324   USE dom_oce        ! ocean space and time domain 
    2425   USE phycst         ! physical constants 
     26   USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files 
     27   ! 
    2528   USE in_out_manager ! I/O manager 
    2629   USE lib_mpp        ! MPP library 
     
    3538 
    3639   !!---------------------------------------------------------------------- 
    37    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     40   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    3841   !! $Id$  
    3942   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    106109      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05 
    107110      INTEGER  ::   isrow                ! index for ORCA1 starting row 
    108  
     111      INTEGER  ::   ie1e2u_v             ! fag for u- & v-surface read in coordinate file or not 
    109112      !!---------------------------------------------------------------------- 
    110113      ! 
     
    122125         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m   
    123126      ENDIF 
    124  
    125  
    126       SELECT CASE( jphgr_msh )   ! type of horizontal mesh 
    127  
    128       CASE ( 0 )                     !  curvilinear coordinate on the sphere read in coordinate.nc file 
    129  
     127      ! 
     128      ! 
     129      SELECT CASE( jphgr_msh )   !  type of horizontal mesh   
     130      ! 
     131      CASE ( 0 )                     !==  read in coordinate.nc file  ==! 
     132         ! 
    130133         IF(lwp) WRITE(numout,*) 
    131134         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file' 
    132  
    133          CALL hgr_read           ! Defaultl option  :   NetCDF file 
    134  
    135          !                                                ! ===================== 
    136          IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    137             !                                             ! ===================== 
    138             IF( nn_cla == 0 ) THEN 
    139                ! 
    140                ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u = 20 km) 
    141                ij0 = 102   ;   ij1 = 102   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3 
    142                IF(lwp) WRITE(numout,*) 
    143                IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar    : e2u reduced to 20 km' 
    144                ! 
    145                ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u = 18 km) 
    146                ij0 =  88   ;   ij1 =  88   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  18.e3 
    147                                                e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  30.e3 
    148                IF(lwp) WRITE(numout,*) 
    149                IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb: e2u reduced to 30 km' 
    150                IF(lwp) WRITE(numout,*) '                                     e1v reduced to 18 km' 
    151             ENDIF 
    152  
    153             ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u = 10 km) 
    154             ij0 = 116   ;   ij1 = 116   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
    155             IF(lwp) WRITE(numout,*) 
    156             IF(lwp) WRITE(numout,*) '             orca_r2: Danish Straits : e2u reduced to 10 km' 
    157             ! 
     135         ! 
     136         ie1e2u_v = 0                  ! set to unread e1e2u and e1e2v 
     137         ! 
     138         CALL hgr_read( ie1e2u_v )     ! read the coordinate.nc file 
     139         ! 
     140         IF( ie1e2u_v == 0 ) THEN      ! e1e2u and e1e2v have not been read: compute them 
     141            !                          ! e2u and e1v does not include a reduction in some strait: apply reduction 
     142            e1e2u (:,:) = e1u(:,:) * e2u(:,:)    
     143            e1e2v (:,:) = e1v(:,:) * e2v(:,:)  
    158144         ENDIF 
    159  
    160             !                                             ! ===================== 
    161          IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
    162             !                                             ! ===================== 
    163             ! This dirty section will be suppressed by simplification process: all this will come back in input files 
    164             ! Currently these hard-wired indices relate to configuration with 
    165             ! extend grid (jpjglo=332) 
    166             ! which had a grid-size of 362x292. 
    167             !  
    168             isrow = 332 - jpjglo 
    169             ! 
    170             ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u = 20 km) 
    171             ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3 
    172             IF(lwp) WRITE(numout,*) 
    173             IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : e2u reduced to 20 km' 
    174  
    175             ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km) 
    176             ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
    177             IF(lwp) WRITE(numout,*) 
    178             IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : e2u reduced to 10 km' 
    179  
    180             ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v = 13 km) 
    181             ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3 
    182             IF(lwp) WRITE(numout,*) 
    183             IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : e1v reduced to 10 km' 
    184  
    185             ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on] 
    186             ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3 
    187             IF(lwp) WRITE(numout,*) 
    188             IF(lwp) WRITE(numout,*) '             orca_r1: Sumba : e1v reduced to 8 km' 
    189  
    190             ii0 =  53           ;   ii1 =  53        ! Ombai Strait (e1v = 13 km) 
    191             ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3 
    192             IF(lwp) WRITE(numout,*) 
    193             IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : e1v reduced to 13 km' 
    194  
    195             ii0 =  56           ;   ii1 =  56        ! Timor Passage (e1v = 20 km) 
    196             ij0 = 164 - isrow   ;   ij1 = 145 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 
    197             IF(lwp) WRITE(numout,*) 
    198             IF(lwp) WRITE(numout,*) '             orca_r1: Timor Passage : e1v reduced to 20 km' 
    199  
    200             ii0 =  55           ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km) 
    201             ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3 
    202             IF(lwp) WRITE(numout,*) 
    203             IF(lwp) WRITE(numout,*) '             orca_r1: W Halmahera : e1v reduced to 30 km' 
    204  
    205             ii0 =  58           ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km) 
    206             ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3 
    207             IF(lwp) WRITE(numout,*) 
    208             IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km' 
    209             ! 
    210             ! 
    211          ENDIF 
    212  
    213          !                                                ! ====================== 
    214          IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration 
    215             !                                             ! ====================== 
    216             ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u = 20 km) 
    217             ij0 = 327   ;   ij1 = 327   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3 
    218             IF(lwp) WRITE(numout,*) 
    219             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Gibraltar Strait' 
    220             ! 
    221             ii0 = 627   ;   ii1 = 628        ! Bosphore Strait (e2u = 10 km) 
    222             ij0 = 343   ;   ij1 = 343   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
    223             IF(lwp) WRITE(numout,*) 
    224             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Bosphore Strait' 
    225             ! 
    226             ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u = 40 km) 
    227             ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  40.e3 
    228             IF(lwp) WRITE(numout,*) 
    229             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Sumba Strait' 
    230             ! 
    231             ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u = 15 km) 
    232             ij0 = 232   ;   ij1 = 232   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  15.e3 
    233             IF(lwp) WRITE(numout,*) 
    234             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Ombai Strait' 
    235             ! 
    236             ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u = 10 km) 
    237             ij0 = 270   ;   ij1 = 270   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
    238             IF(lwp) WRITE(numout,*) 
    239             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e2u at the Palk Strait' 
    240             ! 
    241             ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v = 10 km) 
    242             ij0 = 232   ;   ij1 = 233   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
    243             IF(lwp) WRITE(numout,*) 
    244             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Lombok Strait' 
    245             ! 
    246             ! 
    247             ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v = 25 km) 
    248             ij0 = 276   ;   ij1 = 276   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  25.e3 
    249             IF(lwp) WRITE(numout,*) 
    250             IF(lwp) WRITE(numout,*) '             orca_r05: Reduced e1v at the Bab el Mandeb' 
    251             ! 
    252          ENDIF 
    253  
    254  
    255          ! N.B. :  General case, lat and long function of both i and j indices: 
    256          !     e1t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdila( zti, ztj ) )**2   & 
    257          !                                  + (                           fsdiph( zti, ztj ) )**2  ) 
    258          !     e1u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdila( zui, zuj ) )**2   & 
    259          !                                  + (                           fsdiph( zui, zuj ) )**2  ) 
    260          !     e1v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdila( zvi, zvj ) )**2   & 
    261          !                                  + (                           fsdiph( zvi, zvj ) )**2  ) 
    262          !     e1f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdila( zfi, zfj ) )**2   & 
    263          !                                  + (                           fsdiph( zfi, zfj ) )**2  ) 
    264          ! 
    265          !     e2t(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphit(ji,jj) ) * fsdjla( zti, ztj ) )**2   & 
    266          !                                  + (                           fsdjph( zti, ztj ) )**2  ) 
    267          !     e2u(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiu(ji,jj) ) * fsdjla( zui, zuj ) )**2   & 
    268          !                                  + (                           fsdjph( zui, zuj ) )**2  ) 
    269          !     e2v(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphiv(ji,jj) ) * fsdjla( zvi, zvj ) )**2   & 
    270          !                                  + (                           fsdjph( zvi, zvj ) )**2  ) 
    271          !     e2f(ji,jj) = ra * rad * SQRT(  ( cos( rad*gphif(ji,jj) ) * fsdjla( zfi, zfj ) )**2   & 
    272          !                                  + (                           fsdjph( zfi, zfj ) )**2  ) 
    273  
    274  
    275       CASE ( 1 )                     ! geographical mesh on the sphere with regular grid-spacing 
    276  
     145         ! 
     146      CASE ( 1 )                     !==  geographical mesh on the sphere with regular (in degree) grid-spacing  ==! 
     147         ! 
    277148         IF(lwp) WRITE(numout,*) 
    278149         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing' 
    279150         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg'  
    280  
     151         ! 
    281152         DO jj = 1, jpj 
    282153            DO ji = 1, jpi 
    283                zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - 1 + njmpp - 1 ) 
    284                zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - 1 + njmpp - 1 ) 
    285                zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5 
    286                zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5 
     154               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - 1 + njmpp - 1 ) 
     155               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - 1 + njmpp - 1 ) 
     156               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - 1 + njmpp - 1 ) + 0.5 
     157               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - 1 + njmpp - 1 ) + 0.5 
    287158         ! Longitude 
    288159               glamt(ji,jj) = ppglam0 + ppe1_deg * zti 
     
    307178            END DO 
    308179         END DO 
    309  
    310  
    311       CASE ( 2:3 )                   ! f- or beta-plane with regular grid-spacing 
    312  
     180         ! 
     181      CASE ( 2:3 )                   !==  f- or beta-plane with regular grid-spacing  ==! 
     182         ! 
    313183         IF(lwp) WRITE(numout,*) 
    314184         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing' 
    315185         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'  
    316  
     186         ! 
    317187         ! Position coordinates (in kilometers) 
    318188         !                          ========== 
    319          glam0 = 0.e0 
     189         glam0 = 0._wp 
    320190         gphi0 = - ppe2_m * 1.e-3 
    321           
     191         ! 
    322192#if defined key_agrif  
    323193         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only 
     
    332202         DO jj = 1, jpj 
    333203            DO ji = 1, jpi 
    334                glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 )       ) 
    335                glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( FLOAT( ji - 1 + nimpp - 1 ) + 0.5 ) 
     204               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 )       ) 
     205               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 ) 
    336206               glamv(ji,jj) = glamt(ji,jj) 
    337207               glamf(ji,jj) = glamu(ji,jj) 
    338     
    339                gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 )       ) 
     208               ! 
     209               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 )       ) 
    340210               gphiu(ji,jj) = gphit(ji,jj) 
    341                gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( FLOAT( jj - 1 + njmpp - 1 ) + 0.5 ) 
     211               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 ) 
    342212               gphif(ji,jj) = gphiv(ji,jj) 
    343213            END DO 
    344214         END DO 
    345  
     215         ! 
    346216         ! Horizontal scale factors (in meters) 
    347217         !                              ====== 
     
    350220         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m 
    351221         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m 
    352  
    353       CASE ( 4 )                     ! geographical mesh on the sphere, isotropic MERCATOR type 
    354  
     222         ! 
     223      CASE ( 4 )                     !==  geographical mesh on the sphere, isotropic MERCATOR type  ==! 
     224         ! 
    355225         IF(lwp) WRITE(numout,*) 
    356226         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type' 
    357227         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg' 
    358228         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' ) 
    359  
     229         ! 
    360230         !  Find index corresponding to the equator, given the grid spacing e1_deg 
    361231         !  and the (approximate) southern latitude ppgphi0. 
     
    365235         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg ) 
    366236         IF(  ppgphi0 > 0 )  ijeq = -ijeq 
    367  
     237         ! 
    368238         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq 
    369  
     239         ! 
    370240         DO jj = 1, jpj 
    371241            DO ji = 1, jpi 
    372                zti = FLOAT( ji - 1 + nimpp - 1 )         ;   ztj = FLOAT( jj - ijeq + njmpp - 1 ) 
    373                zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = FLOAT( jj - ijeq + njmpp - 1 ) 
    374                zvi = FLOAT( ji - 1 + nimpp - 1 )         ;   zvj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5 
    375                zfi = FLOAT( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = FLOAT( jj - ijeq + njmpp - 1 ) + 0.5 
     242               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - ijeq + njmpp - 1 ) 
     243               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - ijeq + njmpp - 1 ) 
     244               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 
     245               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 
    376246         ! Longitude 
    377247               glamt(ji,jj) = ppglam0 + ppe1_deg * zti 
     
    396266            END DO 
    397267         END DO 
    398  
    399       CASE ( 5 )                   ! beta-plane with regular grid-spacing and rotated domain (GYRE configuration) 
    400  
     268         ! 
     269      CASE ( 5 )                   !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration) 
     270         ! 
    401271         IF(lwp) WRITE(numout,*) 
    402272         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)' 
    403273         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
    404  
     274         ! 
    405275         ! Position coordinates (in kilometers) 
    406276         !                          ========== 
    407  
     277         ! 
    408278         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN 
    409          zlam1 = -85 
    410          zphi1 = 29 
     279         zlam1 = -85._wp 
     280         zphi1 =  29._wp 
    411281         ! resolution in meters 
    412          ze1 = 106000. / FLOAT(jp_cfg)             
     282         ze1 = 106000. / REAL( jp_cfg , wp )             
    413283         ! benchmark: forced the resolution to be about 100 km 
    414          IF( nbench /= 0 )   ze1 = 106000.e0      
    415          zsin_alpha = - SQRT( 2. ) / 2. 
    416          zcos_alpha =   SQRT( 2. ) / 2. 
     284         IF( nbench /= 0 )   ze1 = 106000._wp      
     285         zsin_alpha = - SQRT( 2._wp ) * 0.5_wp 
     286         zcos_alpha =   SQRT( 2._wp ) * 0.5_wp 
    417287         ze1deg = ze1 / (ra * rad) 
    418          IF( nbench /= 0 )   ze1deg = ze1deg / FLOAT(jp_cfg)        ! benchmark: keep the lat/+lon 
    419          !                                                          ! at the right jp_cfg resolution 
    420          glam0 = zlam1 + zcos_alpha * ze1deg * FLOAT( jpjglo-2 ) 
    421          gphi0 = zphi1 + zsin_alpha * ze1deg * FLOAT( jpjglo-2 ) 
    422  
     288         IF( nbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon 
     289         !                                                           ! at the right jp_cfg resolution 
     290         glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
     291         gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
     292         ! 
    423293         IF( nprint==1 .AND. lwp )   THEN 
    424294            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 
    425295            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0 
    426296         ENDIF 
    427  
     297         ! 
    428298         DO jj = 1, jpj 
    429            DO ji = 1, jpi 
    430              zim1 = FLOAT( ji + nimpp - 1 ) - 1.   ;   zim05 = FLOAT( ji + nimpp - 1 ) - 1.5 
    431              zjm1 = FLOAT( jj + njmpp - 1 ) - 1.   ;   zjm05 = FLOAT( jj + njmpp - 1 ) - 1.5 
    432  
    433              glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
    434              gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
    435  
    436              glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
    437              gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
    438  
    439              glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
    440              gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
    441  
    442              glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
    443              gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
    444            END DO 
    445           END DO 
    446  
     299            DO ji = 1, jpi 
     300               zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5 
     301               zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5 
     302               ! 
     303               glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
     304               gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
     305               ! 
     306               glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
     307               gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
     308               ! 
     309               glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
     310               gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
     311               ! 
     312               glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
     313               gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
     314            END DO 
     315         END DO 
     316         ! 
    447317         ! Horizontal scale factors (in meters) 
    448318         !                              ====== 
     
    451321         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1 
    452322         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1 
    453  
     323         ! 
    454324      CASE DEFAULT 
    455325         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh 
    456326         CALL ctl_stop( ctmp1 ) 
    457  
     327         ! 
    458328      END SELECT 
    459329       
    460       ! T-cell surface 
    461       ! -------------- 
    462       e1e2t(:,:) = e1t(:,:) * e2t(:,:) 
    463      
    464       ! Useful shortcuts (JC: note the duplicated e2e2t array ! Need some cleaning) 
    465       ! --------------------------------------------------------------------------- 
    466       e12t    (:,:) = e1t(:,:) * e2t(:,:) 
    467       e12u    (:,:) = e1u(:,:) * e2u(:,:) 
    468       e12v    (:,:) = e1v(:,:) * e2v(:,:) 
    469       e12f    (:,:) = e1f(:,:) * e2f(:,:) 
    470       r1_e12t (:,:) = 1._wp    / e12t(:,:) 
    471       r1_e12u (:,:) = 1._wp    / e12u(:,:) 
    472       r1_e12v (:,:) = 1._wp    / e12v(:,:) 
    473       r1_e12f (:,:) = 1._wp    / e12f(:,:) 
    474       re2u_e1u(:,:) = e2u(:,:) / e1u(:,:) 
    475       re1v_e2v(:,:) = e1v(:,:) / e2v(:,:) 
    476       r1_e1t  (:,:) = 1._wp    / e1t(:,:) 
    477       r1_e1u  (:,:) = 1._wp    / e1u(:,:) 
    478       r1_e1v  (:,:) = 1._wp    / e1v(:,:) 
    479       r1_e1f  (:,:) = 1._wp    / e1f(:,:) 
    480       r1_e2t  (:,:) = 1._wp    / e2t(:,:) 
    481       r1_e2u  (:,:) = 1._wp    / e2u(:,:) 
    482       r1_e2v  (:,:) = 1._wp    / e2v(:,:) 
    483       r1_e2f  (:,:) = 1._wp    / e2f(:,:) 
    484  
    485       ! Control printing : Grid informations (if not restart) 
    486       ! ---------------- 
    487  
    488       IF( lwp .AND. .NOT.ln_rstart ) THEN 
     330      ! associated horizontal metrics 
     331      ! ----------------------------- 
     332      ! 
     333      r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:) 
     334      r1_e1u(:,:) = 1._wp / e1u(:,:)   ;   r1_e2u (:,:) = 1._wp / e2u(:,:) 
     335      r1_e1v(:,:) = 1._wp / e1v(:,:)   ;   r1_e2v (:,:) = 1._wp / e2v(:,:) 
     336      r1_e1f(:,:) = 1._wp / e1f(:,:)   ;   r1_e2f (:,:) = 1._wp / e2f(:,:) 
     337      ! 
     338      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:) 
     339      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:) 
     340      IF( jphgr_msh /= 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them 
     341         e1e2u (:,:) = e1u(:,:) * e2u(:,:)    
     342         e1e2v (:,:) = e1v(:,:) * e2v(:,:)  
     343      ENDIF 
     344      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in both cases 
     345      r1_e1e2v(:,:) = 1._wp / e1e2v(:,:) 
     346      !    
     347      e2_e1u(:,:) = e2u(:,:) / e1u(:,:) 
     348      e1_e2v(:,:) = e1v(:,:) / e2v(:,:) 
     349 
     350      IF( lwp .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart) 
    489351         WRITE(numout,*) 
    490352         WRITE(numout,*) '          longitude and e1 scale factors' 
     
    4963589300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    & 
    497359            f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 ) 
    498           
     360            ! 
    499361         WRITE(numout,*) 
    500362         WRITE(numout,*) '          latitude and e2 scale factors' 
     
    506368      ENDIF 
    507369 
    508        
    509       IF( nprint == 1 .AND. lwp ) THEN 
    510          WRITE(numout,*) '          e1u e2u ' 
    511          CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    512          CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    513          WRITE(numout,*) '          e1v e2v  ' 
    514          CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    515          CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    516          WRITE(numout,*) '          e1f e2f  ' 
    517          CALL prihre( e1f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    518          CALL prihre( e2f,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    519       ENDIF 
    520  
    521370 
    522371      ! ================= ! 
     
    525374 
    526375      SELECT CASE( jphgr_msh )   ! type of horizontal mesh 
    527  
     376      ! 
    528377      CASE ( 0, 1, 4 )               ! mesh on the sphere 
    529  
     378         ! 
    530379         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) )  
    531  
     380         ! 
    532381      CASE ( 2 )                     ! f-plane at ppgphi0  
    533  
     382         ! 
    534383         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 
    535  
     384         ! 
    536385         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1) 
    537  
     386         ! 
    538387      CASE ( 3 )                     ! beta-plane 
    539  
     388         ! 
    540389         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0 
    541          zphi0   = ppgphi0 - FLOAT( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points 
    542           
     390         zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points 
     391         ! 
    543392#if defined key_agrif 
    544393         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only 
    545394            IF( .NOT. Agrif_Root() ) THEN 
    546               zphi0 = ppgphi0 - FLOAT( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m)   &  
    547                     &           / (ra * rad) 
     395              zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad) 
    548396            ENDIF 
    549397         ENDIF 
    550398#endif          
    551399         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south 
    552  
     400         ! 
    553401         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south) 
    554           
     402         ! 
    555403         IF(lwp) THEN 
    556404            WRITE(numout,*)  
     
    565413            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff 
    566414         END IF 
    567  
     415         ! 
    568416      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration) 
    569  
     417         ! 
    570418         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0 
    571          zphi0 = 15.e0                                                      ! latitude of the first row F-points 
     419         zphi0 = 15._wp                                                     ! latitude of the first row F-points 
    572420         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south 
    573  
     421         ! 
    574422         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south) 
    575  
     423         ! 
    576424         IF(lwp) THEN 
    577425            WRITE(numout,*)  
     
    579427            WRITE(numout,*) '          Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 
    580428         ENDIF 
    581  
     429         ! 
    582430         IF( lk_mpp ) THEN  
    583431            zminff=ff(nldi,nldj) 
     
    587435            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff 
    588436         END IF 
    589  
     437         ! 
    590438      END SELECT 
    591439 
     
    596444 
    597445      IF( nperio == 2 ) THEN 
    598          znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / FLOAT( jpi ) 
     446         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 
    599447         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' ) 
    600448      ENDIF 
     
    605453 
    606454 
    607    SUBROUTINE hgr_read 
     455   SUBROUTINE hgr_read( ke1e2u_v ) 
    608456      !!--------------------------------------------------------------------- 
    609457      !!              ***  ROUTINE hgr_read  *** 
    610458      !! 
    611       !! ** Purpose :   Read a coordinate file in NetCDF format  
    612       !! 
    613       !! ** Method  :   The mesh file has been defined trough a analytical  
    614       !!      or semi-analytical method. It is read in a NetCDF file.  
    615       !!      
     459      !! ** Purpose :   Read a coordinate file in NetCDF format using IOM 
     460      !! 
    616461      !!---------------------------------------------------------------------- 
    617462      USE iom 
    618  
     463      !! 
     464      INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0) 
     465      ! 
    619466      INTEGER ::   inum   ! temporary logical unit 
    620467      !!---------------------------------------------------------------------- 
    621  
     468      ! 
    622469      IF(lwp) THEN 
    623470         WRITE(numout,*) 
     
    625472         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk 
    626473      ENDIF 
    627        
     474      ! 
    628475      CALL iom_open( 'coordinates', inum ) 
    629        
     476      ! 
    630477      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 
    631478      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr ) 
    632479      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr ) 
    633480      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr ) 
    634        
     481      ! 
    635482      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr ) 
    636483      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr ) 
    637484      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr ) 
    638485      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr ) 
    639        
    640       CALL iom_get( inum, jpdom_data, 'e1t', e1t, lrowattr=ln_use_jattr ) 
    641       CALL iom_get( inum, jpdom_data, 'e1u', e1u, lrowattr=ln_use_jattr ) 
    642       CALL iom_get( inum, jpdom_data, 'e1v', e1v, lrowattr=ln_use_jattr ) 
    643       CALL iom_get( inum, jpdom_data, 'e1f', e1f, lrowattr=ln_use_jattr ) 
    644        
    645       CALL iom_get( inum, jpdom_data, 'e2t', e2t, lrowattr=ln_use_jattr ) 
    646       CALL iom_get( inum, jpdom_data, 'e2u', e2u, lrowattr=ln_use_jattr ) 
    647       CALL iom_get( inum, jpdom_data, 'e2v', e2v, lrowattr=ln_use_jattr ) 
    648       CALL iom_get( inum, jpdom_data, 'e2f', e2f, lrowattr=ln_use_jattr ) 
    649        
     486      ! 
     487      CALL iom_get( inum, jpdom_data, 'e1t'  , e1t  , lrowattr=ln_use_jattr ) 
     488      CALL iom_get( inum, jpdom_data, 'e1u'  , e1u  , lrowattr=ln_use_jattr ) 
     489      CALL iom_get( inum, jpdom_data, 'e1v'  , e1v  , lrowattr=ln_use_jattr ) 
     490      CALL iom_get( inum, jpdom_data, 'e1f'  , e1f  , lrowattr=ln_use_jattr ) 
     491      ! 
     492      CALL iom_get( inum, jpdom_data, 'e2t'  , e2t  , lrowattr=ln_use_jattr ) 
     493      CALL iom_get( inum, jpdom_data, 'e2u'  , e2u  , lrowattr=ln_use_jattr ) 
     494      CALL iom_get( inum, jpdom_data, 'e2v'  , e2v  , lrowattr=ln_use_jattr ) 
     495      CALL iom_get( inum, jpdom_data, 'e2f'  , e2f  , lrowattr=ln_use_jattr ) 
     496      ! 
     497      IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN 
     498         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file' 
     499         CALL iom_get( inum, jpdom_data, 'e1e2u'  , e1e2u  , lrowattr=ln_use_jattr ) 
     500         CALL iom_get( inum, jpdom_data, 'e1e2v'  , e1e2v  , lrowattr=ln_use_jattr ) 
     501         ke1e2u_v = 1 
     502      ELSE 
     503         ke1e2u_v = 0 
     504      ENDIF 
     505      ! 
    650506      CALL iom_close( inum ) 
    651507       
     508!!gm   THIS is TO BE REMOVED !!!!!!! 
     509 
    652510! need to be define for the extended grid south of -80S 
    653511! some point are undefined but you need to have e1 and e2 .NE. 0 
     
    676534         e2f=1.0e2 
    677535      END WHERE 
     536!!gm end 
    678537        
    679538    END SUBROUTINE hgr_read 
Note: See TracChangeset for help on using the changeset viewer.