Changeset 6484


Ignore:
Timestamp:
2016-04-19T18:01:48+02:00 (5 years ago)
Author:
flavoni
Message:

update usr_def module, see ticket #1692

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

    r6434 r6484  
    297297         ELSE 
    298298            IF(lwp) WRITE(numout,*) '          ln_read_cfg CALL user_defined module' 
    299             CALL usr_def_hgr( nbench , jp_cfg  ,                    & 
    300             &                 ff,                         & 
    301             &                 glamt  , glamu   , glamv   , glamf   ,         & 
    302             &                 gphit  , gphiu   , gphiv   , gphif   ,         & 
    303             &                 e1t    , e1u     , e1v     , e1f     ,         & 
    304             &                 e2t    , e2u     , e2v     , e2f     ,         &  
    305             &                 e1e2u  , e1e2v   , e2_e1u  , e1_e2v     ) 
     299            CALL usr_def_hgr( nbench , jp_cfg, iff   , ff    ,    & 
     300            &                 glamt  , glamu , glamv , glamf ,    & 
     301            &                 gphit  , gphiu , gphiv , gphif ,    & 
     302            &                 e1t    , e1u   , e1v   , e1f   ,    & 
     303            &                 e2t    , e2u   , e2v   , e2f   ,    &  
     304            &                 ie1e2u_v  ) 
    306305         ! 
    307306         ENDIF 
  • branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90

    r6434 r6484  
    22   !!============================================================================== 
    33   !!                       ***  MODULE usrdef   *** 
    4    !! User defined module: used like example to define init, sbc, bathy, ... 
     4   !! User defined module:  set analytically the minimum information needed to 
     5   !!       define a configuration (domain, initial state, forcing) 
     6   !! 
     7   !!              ===  GYRE configuration given as an example  === 
     8   !!  
     9   !!        This module is given as an example, it can be modified 
     10   !!      following the user's desiderata.  
     11   !!      second order accuracy schemes. 
    512   !!============================================================================== 
    6    !! History :  NEMO ! 2016-03  (S. Flavoni)  
     13   !! History :  NEMO ! 2016-03  (S. Flavoni)  Original code  
    714   !!---------------------------------------------------------------------- 
    815 
    916   !!---------------------------------------------------------------------- 
    1017   !!   usr_def_hgr       : compute the horizontal mesh  
     18   !!   usr_def_zgr       : compute the vertical mesh  
     19   !!  to be defined: 
    1120   !!   usr_def_ini       : initial state  
    1221   !!   usr_def_sbc       : compute the surface bounday conditions 
     
    2130 
    2231   PUBLIC usr_def_hgr 
     32   PUBLIC usr_def_zgr 
    2333 
    2434   !!---------------------------------------------------------------------- 
    25    !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    26    !! $Id: usrdef.F90 6140 2016-03-21 11:35:23Z flavoni $  
     35   !! NEMO/OPA 3.7 , NEMO Consortium (2016) 
     36   !! $Id: $  
    2737   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    2838   !!---------------------------------------------------------------------- 
    2939CONTAINS 
    3040 
    31    SUBROUTINE usr_def_hgr( nbench, jp_cfg, pff,    & 
    32            &              pglamt, pglamu, pglamv , pglamf,     & 
    33            &              pgphit, pgphiu, pgphiv , pgphif,     &           
    34            &              pe1t  , pe1u  , pe1v   , pe1f  ,     &  
    35            &              pe2t  , pe2u  , pe2v   , pe2f  ,     & 
    36            &              pe1e2u, pe1e2v, pe2_e1u, pe1_e2v      ) 
    37       !!---------------------------------------------------------------------- 
     41 
     42   SUBROUTINE usr_def_hgr( nbench, jp_cfg, kff, pff,       & ! Coriolis parameter  (if domain not on the sphere) 
     43           &               pglamt, pglamu, pglamv , pglamf, & ! gridpoints position (required) 
     44           &               pgphit, pgphiu, pgphiv , pgphif, & !          
     45           &               pe1t  , pe1u  , pe1v   , pe1f  , & ! scale factors (required) 
     46           &               pe2t  , pe2u  , pe2v   , pe2f  , & ! 
     47           &               ke1e2u_v                       )! u- & v-surfaces (if gridsize reduction is used in strait(s)) 
     48      !!------------------------------------------------------------------------------------------------- 
    3849      !!                  ***  ROUTINE usr_def_hgr  *** 
    3950      !! 
    40       !! ** Purpose :   compute the geographical position (in degre) of the  
    41       !!      model grid-points, the horizontal scale factors (in meters) and  
    42       !!      the Coriolis factor (in s-1). 
     51      !! ** Purpose :   user defined mesh and Coriolis parameter 
    4352      !! 
    44       !! ** Method  :   the geographical position of the model grid-points is 
    45       !!      defined from analytical functions, fslam and fsphi, the deriva- 
    46       !!      tives of which gives the horizontal scale factors e1,e2. 
    47       !!      Defining two function fslam and fsphi and their derivatives in  
    48       !!      the two horizontal directions (fse1 and fse2), the model grid- 
    49       !!      point position and scale factors are given by: 
    50       !!          t-point: 
    51       !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    ) 
    52       !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    ) 
    53       !!          u-point: 
    54       !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    ) 
    55       !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    ) 
    56       !!          v-point: 
    57       !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2) 
    58       !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2) 
    59       !!          f-point: 
    60       !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2) 
    61       !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2) 
    62       !!      Where fse1 and fse2 are defined by: 
    63       !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2 
    64       !!                                     +          di(fsphi) **2 )(i,j) 
    65       !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2 
    66       !!                                     +          dj(fsphi) **2 )(i,j) 
     53      !! ** Method  :   set all intent(out) argument to a proper value 
    6754      !! 
    68       !!      The coriolis factor is given at z-point by: 
    69       !!         ff = 2.*omega*sin(gphif)      (in s-1) 
     55      !!                Here GYRE configuration : 
     56      !!          Rectangular mid-latitude domain  
     57      !!          - with axes rotated by 45 degrees 
     58      !!          - a constant horizontal resolution of 106 km  
     59      !!          - on a beta-plane 
    7060      !! 
    71       !!      This routine is given as an example, it must be modified 
    72       !!      following the user s desiderata. nevertheless, the output as 
    73       !!      well as the way to compute the model grid-point position and 
    74       !!      horizontal scale factors must be respected in order to insure 
    75       !!      second order accuracy schemes. 
    76       !! 
    77       !! N.B. If the domain is periodic, verify that scale factors are also 
    78       !!      periodic, and the coriolis term again. 
    79       !! 
    80       !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of  
    81       !!                           t-,    u-,   v- and f-points (in degre) 
    82       !!              - define  gphit, gphiu, gphiv, gphif: latitude  of  
    83       !!                           t-,    u-,    v- and f-points (in degre) 
    84       !!              - define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal 
    85       !! scale factors (in meters) at t-,    u-,      v- i and f-points. 
    86       !!              - define ff: coriolis factor at f-point 
    87       !! 
    88       !! References :   Marti, Madec and Delecluse, 1992, JGR 
    89       !!                Madec, Imbard, 1996, Clim. Dyn. 
     61      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees)  
     62      !!              - define coriolis parameter at f-point if the domain in not on the sphere (on beta-plane) 
     63      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters) 
     64      !!              - define u- & v-surfaces (if gridsize reduction is used in some straits) (in m2) 
     65      !!------------------------------------------------------------------------------------------------- 
     66 
    9067      !!---------------------------------------------------------------------- 
    9168      INTEGER                 , INTENT(in   ) ::   nbench, jp_cfg   ! parameter of namelist for benchmark, and dimension of GYRE 
    92       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pff              !  coriolis factor at f-point  
    93       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pglamt, pglamu, pglamv, pglamf !  longitude outputs  
    94       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgphit, pgphiu, pgphiv, pgphif !  latitude outputs 
    95       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe1t, pe1u, pe1v, pe1f   !  horizontal scale factors  
    96       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe2t, pe2u, pe2v, pe2f   !  horizontal scale factors 
    97       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe1e2u , pe1e2v    !  horizontal scale factors 
    98       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe2_e1u, pe1_e2v   !  horizontal scale factors 
    99        !!---------------------------------------------------------------------- 
     69      INTEGER                 , INTENT(  out) ::   kff              ! =1 Coriolis parameter computed here, =0 otherwise 
     70      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pff              ! Coriolis factor at f-point                  [1/s] 
     71      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pglamt, pglamu, pglamv, pglamf !  longitude outputs        [degrees] 
     72      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgphit, pgphiu, pgphiv, pgphif !  latitude outputs         [degrees] 
     73      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe1t, pe1u, pe1v, pe1f   !  i-scale factors                      [m] 
     74      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe2t, pe2u, pe2v, pe2f   !  j-scale factors                      [m] 
     75      INTEGER                 , INTENT(  out) ::   ke1e2u_v                 ! =1 u- & v-surfaces read here, =0 otherwise  
     76      !!------------------------------------------------------------------------------- 
    10077      INTEGER  ::   ji, jj               ! dummy loop indices 
    101       INTEGER  ::   ii0, ii1, ij0, ij1, iff   ! temporary integers 
    102       REAL(wp) ::   zti, zui, zvi, zfi   ! local scalars 
    103       REAL(wp) ::   ztj, zuj, zvj, zfj   !   -      - 
    104       REAL(wp) ::   zphi0, zbeta, znorme ! 
    105       REAL(wp) ::   glam0, gphi0 
    106       REAL(wp) ::   zarg, zf0, zminff, zmaxff 
    107       REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg 
    108       REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05 
    109       !!---------------------------------------------------------------------- 
     78      REAL(wp) ::   zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1  , ze1deg, zf0   ! local scalars 
     79      REAL(wp) ::   zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme        ! local scalars 
     80      !!------------------------------------------------------------------------------- 
    11081      ! 
    11182      IF( nn_timing == 1 )  CALL timing_start('usr_def_hgr') 
     
    12899      IF( nbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon 
    129100      !                                                           ! at the right jp_cfg resolution 
    130       glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
    131       gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
     101      zlam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
     102      zphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
    132103      ! 
    133104      IF( nprint==1 .AND. lwp )   THEN 
    134105         WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 
    135          WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0 
     106         WRITE(numout,*) '          ze1deg', ze1deg, 'zlam0', zlam0, 'zphi0', zphi0 
    136107      ENDIF 
    137108      ! 
     
    143114            !glamt(i,j) longitude at T-point 
    144115            !gphit(i,j) latitude at T-point   
    145             pglamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
    146             pgphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
     116            pglamt(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
     117            pgphit(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
    147118            ! 
    148119            !glamu(i,j) longitude at U-point 
    149120            !gphiu(i,j) latitude at U-point 
    150             pglamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
    151             pgphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
     121            pglamu(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha 
     122            pgphiu(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha 
    152123            ! 
    153124            !glamv(i,j) longitude at V-point 
    154125            !gphiv(i,j) latitude at V-point 
    155             pglamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
    156             pgphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
     126            pglamv(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
     127            pgphiv(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
     128            ! 
    157129            !glamf(i,j) longitude at F-point 
    158130            !gphif(i,j) latitude at F-point  
    159             pglamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
    160             pgphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
     131            pglamf(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha 
     132            pgphif(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha 
    161133         END DO 
    162134      END DO 
    163135      ! 
    164       ! Horizontal scale factors (in meters) 
    165       !                              ====== 
     136      !       !== Horizontal scale factors ==! (in meters) 
     137      !                      
     138      !                             ! constant grid spacing 
    166139      pe1t(:,:) =  ze1     ;      pe2t(:,:) = ze1 
    167140      pe1u(:,:) =  ze1     ;      pe2u(:,:) = ze1 
    168141      pe1v(:,:) =  ze1     ;      pe2v(:,:) = ze1 
    169142      pe1f(:,:) =  ze1     ;      pe2f(:,:) = ze1 
     143      !                             ! NO reduction of grid size in some straits  
     144      ke1e2u_v = 0                  !    ==>> u_ & v_surfaces will be computed in dom_ghr routine 
    170145      ! 
    171       pe1e2u (:,:) = pe1u(:,:) * pe2u(:,:)    
    172       pe1e2v (:,:) = pe1v(:,:) * pe2v(:,:)  
    173146      ! 
    174       pe2_e1u(:,:) = pe2u(:,:) / pe1u(:,:) 
    175       pe1_e2v(:,:) = pe1v(:,:) / pe2v(:,:) 
    176147 
    177       IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart) 
    178          WRITE(numout,*) 
    179          WRITE(numout,*) '          longitude and e1 scale factors' 
    180          WRITE(numout,*) '          ------------------------------' 
    181          WRITE(numout,9300) ( ji, pglamt(ji,1), pglamu(ji,1),   & 
    182             pglamv(ji,1), pglamf(ji,1),   & 
    183             pe1t(ji,1), pe1u(ji,1),   & 
    184             pe1v(ji,1), pe1f(ji,1), ji = 1, jpi,10) 
    185 9300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    & 
    186             f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 ) 
    187             ! 
    188          WRITE(numout,*) 
    189          WRITE(numout,*) '          latitude and e2 scale factors' 
    190          WRITE(numout,*) '          -----------------------------' 
    191          WRITE(numout,9300) ( jj, pgphit(1,jj), pgphiu(1,jj),   & 
    192             &                     pgphiv(1,jj), pgphif(1,jj),   & 
    193             &                     pe2t  (1,jj), pe2u  (1,jj),   & 
    194             &                     pe2v  (1,jj), pe2f  (1,jj), jj = 1, jpj, 10 ) 
    195       ENDIF 
    196  
    197  
    198       ! ================= ! 
    199       !  Coriolis factor  ! 
    200       ! ================= ! 
    201       ! beta-plane and rotated domain (gyre configuration) 
     148      !                      !==  Coriolis parameter  ==! 
     149      kff = 1                                                           !  indicate not to compute ff afterward 
    202150      ! 
    203       !SF old ppsphi0: not more necessary?? zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0 
    204       zbeta = 2. * omega * COS( rad * gphi0 ) / ra                       ! beta at latitude gphi0 
    205       zphi0 = 15._wp                                                     ! latitude of the first row F-points 
    206       zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south 
     151      zbeta = 2. * omega * COS( rad * zphi1 ) / ra                      ! beta at latitude zphi1 
     152      !SF we overwrite zphi0 (south point in latitude) used just above to define pgphif (value of zphi0=15.5190567531966) 
     153      !SF for computation of coriolis we keep the parameter of XXXX 
     154      zphi0 = 15._wp                                                     !  latitude of the most southern grid point   
     155      zf0   = 2. * omega * SIN( rad * zphi0 )                            !  compute f0 1st point south 
    207156      ! 
    208       pff(:,:) = ( zf0 + zbeta * ABS( pgphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south) 
    209       iff = 1 
     157      pff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south) 
    210158      ! 
    211       IF(lwp) THEN 
    212          WRITE(numout,*)  
    213          WRITE(numout,*) '          Beta-plane and rotated domain : ' 
    214          WRITE(numout,*) '          Coriolis parameter varies in this processor from ', pff(nldi,nldj),' to ', pff(nldi,nlej) 
    215       ENDIF 
    216       ! 
    217       IF( lk_mpp ) THEN  
    218          zminff=pff(nldi,nldj) 
    219          zmaxff=pff(nldi,nlej) 
    220          CALL mpp_min( zminff )   ! min over the global domain 
    221          CALL mpp_max( zmaxff )   ! max over the global domain 
    222          IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff 
    223       END IF 
     159      IF(lwp) WRITE(numout,*) '                           beta-plane used. beta = ', zbeta, ' 1/(s.m)' 
    224160      ! 
    225161      IF( nn_timing == 1 )  CALL timing_stop('usr_def_hgr') 
    226162      ! 
    227163   END SUBROUTINE usr_def_hgr 
     164   
     165   SUBROUTINE usr_def_zgr()           ! Coriolis parameter  (if domain not on the sphere) 
     166       ! subroutine for vertical grid 
     167   END SUBROUTINE usr_def_zgr 
    228168 
    229169   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.