Changeset 6593 for branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM
- Timestamp:
- 2016-05-20T16:34:59+02:00 (8 years ago)
- 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 28 28 USE par_oce ! ocean space and time domain 29 29 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 31 32 ! 32 33 USE in_out_manager ! I/O manager 33 34 USE lib_mpp ! MPP library 34 35 USE timing ! Timing 35 USE usrdef ! User defined routine36 36 37 37 IMPLICIT NONE … … 57 57 !! the Coriolis factor (in s-1). 58 58 !! 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) 105 79 !!---------------------------------------------------------------------- 106 80 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 117 92 !!---------------------------------------------------------------------- 118 93 ! … … 122 97 WRITE(numout,*) 123 98 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 ==! 138 104 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 ! 313 120 ! associated horizontal metrics 314 121 ! ----------------------------- … … 321 128 e1e2t (:,:) = e1t(:,:) * e2t(:,:) ; r1_e1e2t(:,:) = 1._wp / e1e2t(:,:) 322 129 e1e2f (:,:) = e1f(:,:) * e2f(:,:) ; r1_e1e2f(:,:) = 1._wp / e1e2f(:,:) 323 IF( jphgr_msh /= 0 ) THEN ! e1e2u and e1e2v have not been set: compute them130 IF( ie1e2u_v == 0 ) THEN ! e1e2u and e1e2v have not been set: compute them 324 131 e1e2u (:,:) = e1u(:,:) * e2u(:,:) 325 132 e1e2v (:,:) = e1v(:,:) * e2v(:,:) … … 330 137 e2_e1u(:,:) = e2u(:,:) / e1u(:,:) 331 138 e1_e2v(:,:) = e1v(:,:) / e2v(:,:) 332 333 139 334 140 ! ================= ! 335 141 ! Coriolis factor ! 336 142 ! ================= ! 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 341 144 ! 342 145 ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 343 146 ! 344 CASE ( 2 ) ! f-plane at ppgphi0345 !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-plane351 !352 zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra ! beta at latitude ppgphi0353 zphi0 = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad ) ! latitude of the first row F-points354 !355 #if defined key_agrif356 IF( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only357 IF( .NOT.Agrif_Root() ) THEN358 zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)359 ENDIF360 ENDIF361 #endif362 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south363 !364 ff(:,:) = ( zf0 + zbeta * gphif(:,:) * 1.e+3 ) ! f = f0 +beta* y ( y=0 at south)365 !366 IF(lwp) THEN367 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 ENDIF371 IF( lk_mpp ) THEN372 zminff=ff(nldi,nldj)373 zmaxff=ff(nldi,nlej)374 CALL mpp_min( zminff ) ! min over the global domain375 CALL mpp_max( zmaxff ) ! max over the global domain376 IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff377 END IF378 !379 END SELECT380 381 382 ! Control of domain for symetrical condition383 ! ------------------------------------------384 ! The equator line must be the latitude coordinate axe385 386 IF( nperio == 2 ) THEN387 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' )389 147 ENDIF 390 148 ! … … 403 161 USE iom 404 162 !! 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) 407 165 ! 408 166 INTEGER :: inum ! temporary logical unit … … 415 173 ENDIF 416 174 ! 417 !SF CALL iom_open( 'coordinates', inum )418 175 CALL iom_open( 'mesh_mask', inum ) 419 176 ! -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90
r6583 r6593 78 78 ! resolution in meters 79 79 ze1 = 106000. / REAL( k_cfg , wp ) 80 ! benchmark: forced the resolution to be about 100 km81 IF( kbench /= 0 ) ze1 = 106000._wp82 80 zsin_alpha = - SQRT( 2._wp ) * 0.5_wp 83 81 zcos_alpha = SQRT( 2._wp ) * 0.5_wp 84 82 ze1deg = ze1 / (ra * rad) 85 IF( kbench /= 0 ) ze1deg = ze1deg / REAL( jp_cfg , wp ) ! benchmark: keep the lat/+lon86 ! ! at the right jp_cfg resolution87 83 zlam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp ) 88 84 zphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp ) 89 85 ! 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 90 88 IF( nprint==1 .AND. lwp ) THEN 91 89 WRITE(numout,*) ' ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
Note: See TracChangeset
for help on using the changeset viewer.