Changeset 6484
- Timestamp:
- 2016-04-19T18:01:48+02:00 (7 years ago)
- Location:
- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 2 edited
-
DOM/domhgr.F90 (modified) (1 diff)
-
usrdef.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r6434 r6484 297 297 ELSE 298 298 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 ) 306 305 ! 307 306 ENDIF -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90
r6434 r6484 2 2 !!============================================================================== 3 3 !! *** 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. 5 12 !!============================================================================== 6 !! History : NEMO ! 2016-03 (S. Flavoni) 13 !! History : NEMO ! 2016-03 (S. Flavoni) Original code 7 14 !!---------------------------------------------------------------------- 8 15 9 16 !!---------------------------------------------------------------------- 10 17 !! usr_def_hgr : compute the horizontal mesh 18 !! usr_def_zgr : compute the vertical mesh 19 !! to be defined: 11 20 !! usr_def_ini : initial state 12 21 !! usr_def_sbc : compute the surface bounday conditions … … 21 30 22 31 PUBLIC usr_def_hgr 32 PUBLIC usr_def_zgr 23 33 24 34 !!---------------------------------------------------------------------- 25 !! NEMO/OPA 3.7 , NEMO Consortium (201 4)26 !! $Id: usrdef.F90 6140 2016-03-21 11:35:23Z flavoni$35 !! NEMO/OPA 3.7 , NEMO Consortium (2016) 36 !! $Id: $ 27 37 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 28 38 !!---------------------------------------------------------------------- 29 39 CONTAINS 30 40 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 !!------------------------------------------------------------------------------------------------- 38 49 !! *** ROUTINE usr_def_hgr *** 39 50 !! 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 43 52 !! 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 67 54 !! 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 70 60 !! 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 90 67 !!---------------------------------------------------------------------- 91 68 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-point93 REAL(wp), DIMENSION(:,:), INTENT( out) :: p glamt, pglamu, pglamv, pglamf ! longitude outputs94 REAL(wp), DIMENSION(:,:), INTENT( out) :: pg phit, pgphiu, pgphiv, pgphif ! latitude outputs95 REAL(wp), DIMENSION(:,:), INTENT( out) :: p e1t, pe1u, pe1v, pe1f ! horizontal scale factors96 REAL(wp), DIMENSION(:,:), INTENT( out) :: pe 2t, pe2u, pe2v, pe2f ! horizontal scale factors97 REAL(wp), DIMENSION(:,:), INTENT( out) :: pe 1e2u , pe1e2v ! horizontal scale factors98 REAL(wp), DIMENSION(:,:), INTENT( out) :: pe2_e1u, pe1_e2v ! horizontal scale factors99 !!----------------------------------------------------------------------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 !!------------------------------------------------------------------------------- 100 77 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 !!------------------------------------------------------------------------------- 110 81 ! 111 82 IF( nn_timing == 1 ) CALL timing_start('usr_def_hgr') … … 128 99 IF( nbench /= 0 ) ze1deg = ze1deg / REAL( jp_cfg , wp ) ! benchmark: keep the lat/+lon 129 100 ! ! 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 ) 132 103 ! 133 104 IF( nprint==1 .AND. lwp ) THEN 134 105 WRITE(numout,*) ' ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 135 WRITE(numout,*) ' ze1deg', ze1deg, ' glam0', glam0, 'gphi0', gphi0106 WRITE(numout,*) ' ze1deg', ze1deg, 'zlam0', zlam0, 'zphi0', zphi0 136 107 ENDIF 137 108 ! … … 143 114 !glamt(i,j) longitude at T-point 144 115 !gphit(i,j) latitude at T-point 145 pglamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha146 pgphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha116 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 147 118 ! 148 119 !glamu(i,j) longitude at U-point 149 120 !gphiu(i,j) latitude at U-point 150 pglamu(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha151 pgphiu(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha121 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 152 123 ! 153 124 !glamv(i,j) longitude at V-point 154 125 !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 ! 157 129 !glamf(i,j) longitude at F-point 158 130 !gphif(i,j) latitude at F-point 159 pglamf(ji,jj) = glam0 + zim1 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha160 pgphif(ji,jj) = gphi0 - zim1 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha131 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 161 133 END DO 162 134 END DO 163 135 ! 164 ! Horizontal scale factors (in meters) 165 ! ====== 136 ! !== Horizontal scale factors ==! (in meters) 137 ! 138 ! ! constant grid spacing 166 139 pe1t(:,:) = ze1 ; pe2t(:,:) = ze1 167 140 pe1u(:,:) = ze1 ; pe2u(:,:) = ze1 168 141 pe1v(:,:) = ze1 ; pe2v(:,:) = ze1 169 142 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 170 145 ! 171 pe1e2u (:,:) = pe1u(:,:) * pe2u(:,:)172 pe1e2v (:,:) = pe1v(:,:) * pe2v(:,:)173 146 ! 174 pe2_e1u(:,:) = pe2u(:,:) / pe1u(:,:)175 pe1_e2v(:,:) = pe1v(:,:) / pe2v(:,:)176 147 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 202 150 ! 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 207 156 ! 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) 210 158 ! 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)' 224 160 ! 225 161 IF( nn_timing == 1 ) CALL timing_stop('usr_def_hgr') 226 162 ! 227 163 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 228 168 229 169 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
