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 6593 for branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM – NEMO

Ignore:
Timestamp:
2016-05-20T16:34:59+02:00 (8 years ago)
Author:
flavoni
Message:

#1692: simplified domhgr and usrdef modules

Location:
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r6579 r6593  
    2828   USE par_oce        ! ocean space and time domain 
    2929   USE phycst         ! physical constants 
    30    USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files 
     30   USE usrdef         ! User defined routine 
     31!   USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files 
    3132   ! 
    3233   USE in_out_manager ! I/O manager 
    3334   USE lib_mpp        ! MPP library 
    3435   USE timing         ! Timing 
    35    USE usrdef        ! User defined routine 
    3636 
    3737   IMPLICIT NONE 
     
    5757      !!      the Coriolis factor (in s-1). 
    5858      !! 
    59       !! ** Method  :   The geographical position of the model grid-points is 
    60       !!      defined from analytical functions, fslam and fsphi, the deriva- 
    61       !!      tives of which gives the horizontal scale factors e1,e2. 
    62       !!      Defining two function fslam and fsphi and their derivatives in  
    63       !!      the two horizontal directions (fse1 and fse2), the model grid- 
    64       !!      point position and scale factors are given by: 
    65       !!         t-point: 
    66       !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    ) 
    67       !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    ) 
    68       !!         u-point: 
    69       !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    ) 
    70       !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    ) 
    71       !!         v-point: 
    72       !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2) 
    73       !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2) 
    74       !!            f-point: 
    75       !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2) 
    76       !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2) 
    77       !!      Where fse1 and fse2 are defined by: 
    78       !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2 
    79       !!                                     +          di(fsphi) **2 )(i,j) 
    80       !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2 
    81       !!                                     +          dj(fsphi) **2 )(i,j) 
    82       !! 
    83       !!        The coriolis factor is given at z-point by: 
    84       !!                     ff = 2.*omega*sin(gphif)      (in s-1) 
    85       !! 
    86       !!        This routine is given as an example, it must be modified 
    87       !!      following the user s desiderata. nevertheless, the output as 
    88       !!      well as the way to compute the model grid-point position and 
    89       !!      horizontal scale factors must be respected in order to insure 
    90       !!      second order accuracy schemes. 
    91       !! 
    92       !! N.B. If the domain is periodic, verify that scale factors are also 
    93       !!      periodic, and the coriolis term again. 
    94       !! 
    95       !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-,  
    96       !!                u-, v- and f-points (in degre) 
    97       !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-, 
    98       !!               u-, v-  and f-points (in degre) 
    99       !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal 
    100       !!      scale factors (in meters) at t-, u-, v-, and f-points. 
    101       !!        define ff: coriolis factor at f-point 
    102       !! 
    103       !! References :   Marti, Madec and Delecluse, 1992, JGR 
    104       !!                Madec, Imbard, 1996, Clim. Dyn. 
     59      !! ** Method  :   Controlled by ln_read_cfg logical 
     60      !!              =T : all needed arrays are read in mesh_mask.nc file  
     61      !!              =F : user-defined configuration, all needed arrays  
     62      !!                   are computed in usr-def_hgr subroutine  
     63      !! 
     64      !!                If Coriolis factor is neither read nor computed (iff=0) 
     65      !!              it is computed from gphit assuming that the mesh is 
     66      !!              defined on the sphere : 
     67      !!                   ff = 2.*omega*sin(gphif)      (in s-1) 
     68      !!     
     69      !!                If u- & v-surfaces are neither read nor computed (ie1e2u_v=0) 
     70      !!              (i.e. no use of reduced scale factors in some straits) 
     71      !!              they are computed from e1u, e2u, e1v and e2v as: 
     72      !!                   e1e2u = e1u*e2u   and   e1e2v = e1v*e2v   
     73      !!     
     74      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees) 
     75      !!              - define Coriolis parameter at f-point                   (in 1/s) 
     76      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters) 
     77      !!              - define associated horizontal metrics at t-, u-, v- and f-points 
     78      !!                (inverse of scale factors 1/e1 & 1/e2, surface e1*e2, ratios e1/e2 & e2/e1) 
    10579      !!---------------------------------------------------------------------- 
    10680      INTEGER  ::   ji, jj                   ! dummy loop indices 
    107       INTEGER  ::   ii0, ii1, ij0, ij1, iff  ! temporary integers 
    108       INTEGER  ::   ijeq                     ! index of equator T point (used in case 4) 
    109       REAL(wp) ::   zti, zui, zvi, zfi       ! local scalars 
    110       REAL(wp) ::   ztj, zuj, zvj, zfj       !   -      - 
    111       REAL(wp) ::   zphi0, zbeta, znorme     ! 
    112       REAL(wp) ::   zarg, zf0, zminff, zmaxff 
    113       REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg 
    114       REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05 
    115       INTEGER  ::   isrow                    ! index for ORCA1 starting row 
    116       INTEGER  ::   ie1e2u_v                 ! fag for u- & v-surface read in coordinate file or not 
     81      !INTEGER  ::   ii0, ii1, ij0, ij1, iff  ! temporary integers 
     82      !INTEGER  ::   ijeq                     ! index of equator T point (used in case 4) 
     83      !REAL(wp) ::   zti, zui, zvi, zfi       ! local scalars 
     84      !REAL(wp) ::   ztj, zuj, zvj, zfj       !   -      - 
     85      !REAL(wp) ::   zphi0, zbeta, znorme     ! 
     86      !REAL(wp) ::   zarg, zf0, zminff, zmaxff 
     87      !REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg 
     88      !REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05 
     89      !INTEGER  ::   isrow                    ! index for ORCA1 starting row 
     90      INTEGER  ::   ie1e2u_v                 ! flag for u- & v-surfaces  
     91      INTEGER  ::   iff                      ! flag for Coriolis parameter 
    11792      !!---------------------------------------------------------------------- 
    11893      ! 
     
    12297         WRITE(numout,*) 
    12398         WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters ' 
    124          WRITE(numout,*) '~~~~~~~      type of horizontal mesh           jphgr_msh = ', jphgr_msh 
    125          WRITE(numout,*) '             position of the first row and     ppglam0  = ', ppglam0 
    126          WRITE(numout,*) '             column grid-point (degrees)       ppgphi0  = ', ppgphi0 
    127          WRITE(numout,*) '             zonal      grid-spacing (degrees) ppe1_deg = ', ppe1_deg 
    128          WRITE(numout,*) '             meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg 
    129          WRITE(numout,*) '             zonal      grid-spacing (meters)  ppe1_m   = ', ppe1_m   
    130          WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m   
    131       ENDIF 
    132       ! 
    133       ! 
    134       SELECT CASE( jphgr_msh )   !  type of horizontal mesh   
    135       ! 
    136       CASE ( 0 )                     !==  read in coordinate.nc file  ==! 
    137          ! 
     99         WRITE(numout,*) '~~~~~~~   ' 
     100      ENDIF 
     101      ! 
     102      ! 
     103      IF( ln_read_cfg ) THEN        !==  read in mesh_mask.nc file  ==! 
    138104         IF(lwp) WRITE(numout,*) 
    139          IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file' 
    140          ! 
    141          ie1e2u_v = 0             ! set to unread e1e2u and e1e2v 
    142          iff = 0                  ! set to unread iff 
    143          ! 
    144          CALL hgr_read( ie1e2u_v, iff )     ! read the coordinate.nc file !SF to be changed in mesh_mask 
    145          ! 
    146          IF( ie1e2u_v == 0 ) THEN      ! e1e2u and e1e2v have not been read: compute them 
    147             !                          ! e2u and e1v does not include a reduction in some strait: apply reduction 
    148             e1e2u (:,:) = e1u(:,:) * e2u(:,:)    
    149             e1e2v (:,:) = e1v(:,:) * e2v(:,:)  
    150          ENDIF 
    151          ! 
    152       CASE ( 1 )                     !==  geographical mesh on the sphere with regular (in degree) grid-spacing  ==! 
    153          ! 
    154          IF(lwp) WRITE(numout,*) 
    155          IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing' 
    156          IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg'  
    157          ! 
    158          DO jj = 1, jpj 
    159             DO ji = 1, jpi 
    160                zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - 1 + njmpp - 1 ) 
    161                zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - 1 + njmpp - 1 ) 
    162                zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - 1 + njmpp - 1 ) + 0.5 
    163                zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - 1 + njmpp - 1 ) + 0.5 
    164          ! Longitude 
    165                glamt(ji,jj) = ppglam0 + ppe1_deg * zti 
    166                glamu(ji,jj) = ppglam0 + ppe1_deg * zui 
    167                glamv(ji,jj) = ppglam0 + ppe1_deg * zvi 
    168                glamf(ji,jj) = ppglam0 + ppe1_deg * zfi 
    169          ! Latitude 
    170                gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj 
    171                gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj 
    172                gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj 
    173                gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj 
    174          ! e1 
    175                e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg 
    176                e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg 
    177                e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg 
    178                e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg 
    179          ! e2 
    180                e2t(ji,jj) = ra * rad * ppe2_deg 
    181                e2u(ji,jj) = ra * rad * ppe2_deg 
    182                e2v(ji,jj) = ra * rad * ppe2_deg 
    183                e2f(ji,jj) = ra * rad * ppe2_deg 
    184             END DO 
    185          END DO 
    186          ! 
    187       CASE ( 2:3 )                   !==  f- or beta-plane with regular grid-spacing  ==! 
    188          ! 
    189          IF(lwp) WRITE(numout,*) 
    190          IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing' 
    191          IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'  
    192          ! 
    193          ! Position coordinates (in kilometers) 
    194          !                          ========== 
    195          glam0 = 0._wp 
    196          gphi0 = - ppe2_m * 1.e-3 
    197          ! 
    198 #if defined key_agrif  
    199          IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only 
    200             IF( .NOT. Agrif_Root() ) THEN 
    201               glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3 
    202               gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3 
    203               ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox() 
    204               ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()           
    205             ENDIF 
    206          ENDIF 
    207 #endif          
    208          DO jj = 1, jpj 
    209             DO ji = 1, jpi 
    210                glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 )       ) 
    211                glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 ) 
    212                glamv(ji,jj) = glamt(ji,jj) 
    213                glamf(ji,jj) = glamu(ji,jj) 
    214                ! 
    215                gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 )       ) 
    216                gphiu(ji,jj) = gphit(ji,jj) 
    217                gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 ) 
    218                gphif(ji,jj) = gphiv(ji,jj) 
    219             END DO 
    220          END DO 
    221          ! 
    222          ! Horizontal scale factors (in meters) 
    223          !                              ====== 
    224          e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m 
    225          e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m 
    226          e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m 
    227          e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m 
    228          ! 
    229       CASE ( 4 )                     !==  geographical mesh on the sphere, isotropic MERCATOR type  ==! 
    230          ! 
    231          IF(lwp) WRITE(numout,*) 
    232          IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type' 
    233          IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg' 
    234          IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' ) 
    235          ! 
    236          !  Find index corresponding to the equator, given the grid spacing e1_deg 
    237          !  and the (approximate) southern latitude ppgphi0. 
    238          !  This way we ensure that the equator is at a "T / U" point, when in the domain. 
    239          !  The formula should work even if the equator is outside the domain. 
    240          zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2. 
    241          ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg ) 
    242          IF(  ppgphi0 > 0 )  ijeq = -ijeq 
    243          ! 
    244          IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq 
    245          ! 
    246          DO jj = 1, jpj 
    247             DO ji = 1, jpi 
    248                zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - ijeq + njmpp - 1 ) 
    249                zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - ijeq + njmpp - 1 ) 
    250                zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 
    251                zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5 
    252          ! Longitude 
    253                glamt(ji,jj) = ppglam0 + ppe1_deg * zti 
    254                glamu(ji,jj) = ppglam0 + ppe1_deg * zui 
    255                glamv(ji,jj) = ppglam0 + ppe1_deg * zvi 
    256                glamf(ji,jj) = ppglam0 + ppe1_deg * zfi 
    257          ! Latitude 
    258                gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) ) 
    259                gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) ) 
    260                gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) ) 
    261                gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) ) 
    262          ! e1 
    263                e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg 
    264                e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg 
    265                e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg 
    266                e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg 
    267          ! e2 
    268                e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg 
    269                e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg 
    270                e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg 
    271                e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg 
    272             END DO 
    273          END DO 
    274          ! 
    275       CASE ( 5 )                   !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration) 
    276          ! 
    277          IF(lwp) WRITE(numout,*) 
    278          IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)' 
    279          IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
    280          ! 
    281          IF(lwp) THEN 
    282             WRITE(numout,*) '               control print of value of ln_read_cfg     = ', ln_read_cfg 
    283          ENDIF 
    284          IF (ln_read_cfg) THEN  
    285             IF(lwp) WRITE(numout,*) '          ln_read_cfg read input files' 
    286             ie1e2u_v = 0                  ! set to unread e1e2u and e1e2v 
    287             iff = 0                       ! set to unread ff 
    288             ! 
    289             CALL hgr_read( ie1e2u_v, iff )     ! read the coordinate.nc file !SF to be changed in mesh_mask 
    290             ! 
    291             IF( ie1e2u_v == 0 ) THEN   ! e1e2u and e1e2v have not been read: compute them 
    292                !                       ! e2u and e1v does not include a reduction in some strait: apply reduction 
    293                e1e2u (:,:) = e1u(:,:) * e2u(:,:) 
    294                e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
    295             ENDIF 
    296          ELSE 
    297             IF(lwp) WRITE(numout,*) '          ln_read_cfg CALL user_defined module' 
    298             CALL usr_def_hgr( nbench , jp_cfg, iff   , ff    ,    & 
    299             &                 glamt  , glamu , glamv , glamf ,    & 
    300             &                 gphit  , gphiu , gphiv , gphif ,    & 
    301             &                 e1t    , e1u   , e1v   , e1f   ,    & 
    302             &                 e2t    , e2u   , e2v   , e2f   ,    & 
    303             &                 ie1e2u_v  ) 
    304          ! 
    305          ENDIF 
    306          ! 
    307       CASE DEFAULT 
    308          WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh 
    309          CALL ctl_stop( ctmp1 ) 
    310          ! 
    311       END SELECT 
    312        
     105         IF(lwp) WRITE(numout,*) '          read horizontal mesh in "mesh_mask" file' 
     106         ! 
     107         CALL hgr_read( ie1e2u_v, iff )      
     108         ! 
     109      ELSE     
     110         IF(lwp) WRITE(numout,*) '          ln_read_cfg CALL user_defined module' 
     111         CALL usr_def_hgr( nbench , jp_cfg, iff   , ff    ,    & 
     112         &                 glamt  , glamu , glamv , glamf ,    & 
     113         &                 gphit  , gphiu , gphiv , gphif ,    & 
     114         &                 e1t    , e1u   , e1v   , e1f   ,    & 
     115         &                 e2t    , e2u   , e2v   , e2f   ,    & 
     116         &                 ie1e2u_v  ) 
     117         ! 
     118      ENDIF 
     119      ! 
    313120      ! associated horizontal metrics 
    314121      ! ----------------------------- 
     
    321128      e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:) 
    322129      e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:) 
    323       IF( jphgr_msh /= 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them 
     130      IF( ie1e2u_v == 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them 
    324131         e1e2u (:,:) = e1u(:,:) * e2u(:,:)    
    325132         e1e2v (:,:) = e1v(:,:) * e2v(:,:)  
     
    330137      e2_e1u(:,:) = e2u(:,:) / e1u(:,:) 
    331138      e1_e2v(:,:) = e1v(:,:) / e2v(:,:) 
    332  
    333  
     139       
    334140      ! ================= ! 
    335141      !  Coriolis factor  ! 
    336142      ! ================= ! 
    337  
    338       SELECT CASE( jphgr_msh )   ! type of horizontal mesh 
    339       ! 
    340       CASE ( 0, 1, 4 )               ! mesh on the sphere 
     143      IF ( iff == 0 ) THEN       ! Coriolis parameter has not been defined: compute it on the sphere  
    341144         ! 
    342145         ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) )  
    343146         ! 
    344       CASE ( 2 )                     ! f-plane at ppgphi0  
    345          ! 
    346          ff(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 
    347          ! 
    348          IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1) 
    349          ! 
    350       CASE ( 3 )                     ! beta-plane 
    351          ! 
    352          zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0 
    353          zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points 
    354          ! 
    355 #if defined key_agrif 
    356          IF( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN       ! for EEL6 configuration only 
    357             IF( .NOT.Agrif_Root() ) THEN 
    358               zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad) 
    359             ENDIF 
    360          ENDIF 
    361 #endif          
    362          zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south 
    363          ! 
    364          ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south) 
    365          ! 
    366          IF(lwp) THEN 
    367             WRITE(numout,*)  
    368             WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj) 
    369             WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 
    370          ENDIF 
    371          IF( lk_mpp ) THEN  
    372             zminff=ff(nldi,nldj) 
    373             zmaxff=ff(nldi,nlej) 
    374             CALL mpp_min( zminff )   ! min over the global domain 
    375             CALL mpp_max( zmaxff )   ! max over the global domain 
    376             IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff 
    377          END IF 
    378          ! 
    379       END SELECT 
    380  
    381  
    382       ! Control of domain for symetrical condition 
    383       ! ------------------------------------------ 
    384       ! The equator line must be the latitude coordinate axe 
    385  
    386       IF( nperio == 2 ) THEN 
    387          znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 
    388          IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' ) 
    389147      ENDIF 
    390148      ! 
     
    403161      USE iom 
    404162      !! 
    405       INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in mesh_mask file (=1) or not (=0) 
    406       INTEGER, INTENT( inout ) ::   kff        ! fag: kff read in mesh_mask file (=1) or not (=0) 
     163      INTEGER, INTENT( out ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in mesh_mask file (=1) or not (=0) 
     164      INTEGER, INTENT( out ) ::   kff        ! fag: kff read in mesh_mask file (=1) or not (=0) 
    407165      ! 
    408166      INTEGER ::   inum   ! temporary logical unit 
     
    415173      ENDIF 
    416174      ! 
    417       !SF CALL iom_open( 'coordinates', inum ) 
    418175      CALL iom_open( 'mesh_mask', inum ) 
    419176      ! 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90

    r6583 r6593  
    7878      ! resolution in meters 
    7979      ze1 = 106000. / REAL( k_cfg , wp ) 
    80       ! benchmark: forced the resolution to be about 100 km 
    81             IF( kbench /= 0 )   ze1 = 106000._wp    
    8280      zsin_alpha = - SQRT( 2._wp ) * 0.5_wp 
    8381      zcos_alpha =   SQRT( 2._wp ) * 0.5_wp 
    8482      ze1deg = ze1 / (ra * rad) 
    85       IF( kbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon 
    86       !                                                           ! at the right jp_cfg resolution 
    8783      zlam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
    8884      zphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
    8985      !    
     86      IF( kbench == 1 )  ze1 = 106000._wp ! benchmark: forced the resolution to be about 100 km   
     87      !                                   ! but keep (lat,lon) at the right jp_cfg resolution 
    9088      IF( nprint==1 .AND. lwp )   THEN 
    9189         WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 
Note: See TracChangeset for help on using the changeset viewer.