Changeset 7646 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r6140 r7646 16 16 !! 3.7 ! 2015-09 (G. Madec, S. Flavoni) add cell surface and their inverse 17 17 !! add optional read of e1e2u & e1e2v 18 !! - ! 2016-04 (S. Flavoni, G. Madec) new configuration interface: read or usrdef.F90 18 19 !!---------------------------------------------------------------------- 19 20 20 21 !!---------------------------------------------------------------------- 21 22 !! dom_hgr : initialize the horizontal mesh 22 !! hgr_read : read "coordinate" NetCDFfile23 !! hgr_read : read horizontal information in the domain configuration file 23 24 !!---------------------------------------------------------------------- 24 25 USE dom_oce ! ocean space and time domain 26 USE par_oce ! ocean space and time domain 25 27 USE phycst ! physical constants 26 USE domwri ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files28 USE usrdef_hgr ! User defined routine 27 29 ! 28 30 USE in_out_manager ! I/O manager 31 USE iom ! I/O library 29 32 USE lib_mpp ! MPP library 30 33 USE timing ! Timing … … 33 36 PRIVATE 34 37 35 REAL(wp) :: glam0, gphi0 ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce36 37 38 PUBLIC dom_hgr ! called by domain.F90 38 39 39 40 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 3.7 , NEMO Consortium (201 4)41 !! NEMO/OPA 3.7 , NEMO Consortium (2016) 41 42 !! $Id$ 42 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 48 49 !! *** ROUTINE dom_hgr *** 49 50 !! 50 !! ** Purpose : Compute the geographical position (in degre) of the 51 !! model grid-points, the horizontal scale factors (in meters) and 52 !! the Coriolis factor (in s-1). 53 !! 54 !! ** Method : The geographical position of the model grid-points is 55 !! defined from analytical functions, fslam and fsphi, the deriva- 56 !! tives of which gives the horizontal scale factors e1,e2. 57 !! Defining two function fslam and fsphi and their derivatives in 58 !! the two horizontal directions (fse1 and fse2), the model grid- 59 !! point position and scale factors are given by: 60 !! t-point: 61 !! glamt(i,j) = fslam(i ,j ) e1t(i,j) = fse1(i ,j ) 62 !! gphit(i,j) = fsphi(i ,j ) e2t(i,j) = fse2(i ,j ) 63 !! u-point: 64 !! glamu(i,j) = fslam(i+1/2,j ) e1u(i,j) = fse1(i+1/2,j ) 65 !! gphiu(i,j) = fsphi(i+1/2,j ) e2u(i,j) = fse2(i+1/2,j ) 66 !! v-point: 67 !! glamv(i,j) = fslam(i ,j+1/2) e1v(i,j) = fse1(i ,j+1/2) 68 !! gphiv(i,j) = fsphi(i ,j+1/2) e2v(i,j) = fse2(i ,j+1/2) 69 !! f-point: 70 !! glamf(i,j) = fslam(i+1/2,j+1/2) e1f(i,j) = fse1(i+1/2,j+1/2) 71 !! gphif(i,j) = fsphi(i+1/2,j+1/2) e2f(i,j) = fse2(i+1/2,j+1/2) 72 !! Where fse1 and fse2 are defined by: 73 !! fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2 74 !! + di(fsphi) **2 )(i,j) 75 !! fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2 76 !! + dj(fsphi) **2 )(i,j) 77 !! 78 !! The coriolis factor is given at z-point by: 79 !! ff = 2.*omega*sin(gphif) (in s-1) 80 !! 81 !! This routine is given as an example, it must be modified 82 !! following the user s desiderata. nevertheless, the output as 83 !! well as the way to compute the model grid-point position and 84 !! horizontal scale factors must be respected in order to insure 85 !! second order accuracy schemes. 86 !! 87 !! N.B. If the domain is periodic, verify that scale factors are also 88 !! periodic, and the coriolis term again. 89 !! 90 !! ** Action : - define glamt, glamu, glamv, glamf: longitude of t-, 91 !! u-, v- and f-points (in degre) 92 !! - define gphit, gphiu, gphiv, gphit: latitude of t-, 93 !! u-, v- and f-points (in degre) 94 !! define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal 95 !! scale factors (in meters) at t-, u-, v-, and f-points. 96 !! define ff: coriolis factor at f-point 97 !! 98 !! References : Marti, Madec and Delecluse, 1992, JGR 99 !! Madec, Imbard, 1996, Clim. Dyn. 100 !!---------------------------------------------------------------------- 101 INTEGER :: ji, jj ! dummy loop indices 102 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers 103 INTEGER :: ijeq ! index of equator T point (used in case 4) 104 REAL(wp) :: zti, zui, zvi, zfi ! local scalars 105 REAL(wp) :: ztj, zuj, zvj, zfj ! - - 106 REAL(wp) :: zphi0, zbeta, znorme ! 107 REAL(wp) :: zarg, zf0, zminff, zmaxff 108 REAL(wp) :: zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg 109 REAL(wp) :: zphi1, zsin_alpha, zim05, zjm05 110 INTEGER :: isrow ! index for ORCA1 starting row 111 INTEGER :: ie1e2u_v ! fag for u- & v-surface read in coordinate file or not 51 !! ** Purpose : Read or compute the geographical position (in degrees) 52 !! of the model grid-points, the horizontal scale factors (in meters), 53 !! the associated horizontal metrics, and the Coriolis factor (in s-1). 54 !! 55 !! ** Method : Controlled by ln_read_cfg logical 56 !! =T : all needed arrays are read in mesh_mask.nc file 57 !! =F : user-defined configuration, all needed arrays 58 !! are computed in usr-def_hgr subroutine 59 !! 60 !! If Coriolis factor is neither read nor computed (iff=0) 61 !! it is computed from gphit assuming that the mesh is 62 !! defined on the sphere : 63 !! ff = 2.*omega*sin(gphif) (in s-1) 64 !! 65 !! If u- & v-surfaces are neither read nor computed (ie1e2u_v=0) 66 !! (i.e. no use of reduced scale factors in some straits) 67 !! they are computed from e1u, e2u, e1v and e2v as: 68 !! e1e2u = e1u*e2u and e1e2v = e1v*e2v 69 !! 70 !! ** Action : - define longitude & latitude of t-, u-, v- and f-points (in degrees) 71 !! - define Coriolis parameter at f-point (in 1/s) 72 !! - define i- & j-scale factors at t-, u-, v- and f-points (in meters) 73 !! - define associated horizontal metrics at t-, u-, v- and f-points 74 !! (inverse of scale factors 1/e1 & 1/e2, surface e1*e2, ratios e1/e2 & e2/e1) 75 !!---------------------------------------------------------------------- 76 INTEGER :: ji, jj ! dummy loop indices 77 INTEGER :: ie1e2u_v ! flag for u- & v-surfaces 78 INTEGER :: iff ! flag for Coriolis parameter 112 79 !!---------------------------------------------------------------------- 113 80 ! … … 117 84 WRITE(numout,*) 118 85 WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters ' 119 WRITE(numout,*) '~~~~~~~ type of horizontal mesh jphgr_msh = ', jphgr_msh 120 WRITE(numout,*) ' position of the first row and ppglam0 = ', ppglam0 121 WRITE(numout,*) ' column grid-point (degrees) ppgphi0 = ', ppgphi0 122 WRITE(numout,*) ' zonal grid-spacing (degrees) ppe1_deg = ', ppe1_deg 123 WRITE(numout,*) ' meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg 124 WRITE(numout,*) ' zonal grid-spacing (meters) ppe1_m = ', ppe1_m 125 WRITE(numout,*) ' meridional grid-spacing (meters) ppe2_m = ', ppe2_m 126 ENDIF 127 ! 128 ! 129 SELECT CASE( jphgr_msh ) ! type of horizontal mesh 130 ! 131 CASE ( 0 ) !== read in coordinate.nc file ==! 132 ! 86 WRITE(numout,*) '~~~~~~~ ' 87 WRITE(numout,*) ' namcfg : read (=T) or user defined (=F) configuration ln_read_cfg = ', ln_read_cfg 88 ENDIF 89 ! 90 ! 91 IF( ln_read_cfg ) THEN !== read in mesh_mask.nc file ==! 133 92 IF(lwp) WRITE(numout,*) 134 IF(lwp) WRITE(numout,*) ' curvilinear coordinate on the sphere read in "coordinate" file' 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(:,:) 93 IF(lwp) WRITE(numout,*) ' read horizontal mesh in ', TRIM( cn_domcfg ), ' file' 94 ! 95 CALL hgr_read ( glamt , glamu , glamv , glamf , & ! geographic position (required) 96 & gphit , gphiu , gphiv , gphif , & ! - - 97 & iff , ff_f , ff_t , & ! Coriolis parameter (if not on the sphere) 98 & e1t , e1u , e1v , e1f , & ! scale factors (required) 99 & e2t , e2u , e2v , e2f , & ! - - - 100 & ie1e2u_v , e1e2u , e1e2v ) ! u- & v-surfaces (if gridsize reduction in some straits) 101 ! 102 ELSE !== User defined configuration ==! 103 IF(lwp) WRITE(numout,*) 104 IF(lwp) WRITE(numout,*) ' User defined horizontal mesh (usr_def_hgr)' 105 ! 106 CALL usr_def_hgr( glamt , glamu , glamv , glamf , & ! geographic position (required) 107 & gphit , gphiu , gphiv , gphif , & ! 108 & iff , ff_f , ff_t , & ! Coriolis parameter (if domain not on the sphere) 109 & e1t , e1u , e1v , e1f , & ! scale factors (required) 110 & e2t , e2u , e2v , e2f , & ! 111 & ie1e2u_v , e1e2u , e1e2v ) ! u- & v-surfaces (if gridsize reduction is used in strait(s)) 112 ! 113 ENDIF 114 ! 115 ! !== Coriolis parameter ==! (if necessary) 116 ! 117 IF( iff == 0 ) THEN ! Coriolis parameter has not been defined 118 IF(lwp) WRITE(numout,*) ' Coriolis parameter calculated on the sphere from gphif & gphit' 119 ff_f(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) ! compute it on the sphere at f-point 120 ff_t(:,:) = 2. * omega * SIN( rad * gphit(:,:) ) ! - - - at t-point 121 ELSE 122 IF( ln_read_cfg ) THEN 123 IF(lwp) WRITE(numout,*) ' Coriolis parameter have been read in ', TRIM( cn_domcfg ), ' file' 124 ELSE 125 IF(lwp) WRITE(numout,*) ' Coriolis parameter have been set in usr_def_hgr routine' 144 126 ENDIF 145 ! 146 CASE ( 1 ) !== geographical mesh on the sphere with regular (in degree) grid-spacing ==! 147 ! 148 IF(lwp) WRITE(numout,*) 149 IF(lwp) WRITE(numout,*) ' geographical mesh on the sphere with regular grid-spacing' 150 IF(lwp) WRITE(numout,*) ' given by ppe1_deg and ppe2_deg' 151 ! 152 DO jj = 1, jpj 153 DO ji = 1, jpi 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 158 ! Longitude 159 glamt(ji,jj) = ppglam0 + ppe1_deg * zti 160 glamu(ji,jj) = ppglam0 + ppe1_deg * zui 161 glamv(ji,jj) = ppglam0 + ppe1_deg * zvi 162 glamf(ji,jj) = ppglam0 + ppe1_deg * zfi 163 ! Latitude 164 gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj 165 gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj 166 gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj 167 gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj 168 ! e1 169 e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg 170 e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg 171 e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg 172 e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg 173 ! e2 174 e2t(ji,jj) = ra * rad * ppe2_deg 175 e2u(ji,jj) = ra * rad * ppe2_deg 176 e2v(ji,jj) = ra * rad * ppe2_deg 177 e2f(ji,jj) = ra * rad * ppe2_deg 178 END DO 179 END DO 180 ! 181 CASE ( 2:3 ) !== f- or beta-plane with regular grid-spacing ==! 182 ! 183 IF(lwp) WRITE(numout,*) 184 IF(lwp) WRITE(numout,*) ' f- or beta-plane with regular grid-spacing' 185 IF(lwp) WRITE(numout,*) ' given by ppe1_m and ppe2_m' 186 ! 187 ! Position coordinates (in kilometers) 188 ! ========== 189 glam0 = 0._wp 190 gphi0 = - ppe2_m * 1.e-3 191 ! 192 #if defined key_agrif 193 IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only 194 IF( .NOT. Agrif_Root() ) THEN 195 glam0 = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3 196 gphi0 = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3 197 ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox() 198 ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy() 199 ENDIF 200 ENDIF 201 #endif 202 DO jj = 1, jpj 203 DO ji = 1, jpi 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 ) 206 glamv(ji,jj) = glamt(ji,jj) 207 glamf(ji,jj) = glamu(ji,jj) 208 ! 209 gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) ) 210 gphiu(ji,jj) = gphit(ji,jj) 211 gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 ) 212 gphif(ji,jj) = gphiv(ji,jj) 213 END DO 214 END DO 215 ! 216 ! Horizontal scale factors (in meters) 217 ! ====== 218 e1t(:,:) = ppe1_m ; e2t(:,:) = ppe2_m 219 e1u(:,:) = ppe1_m ; e2u(:,:) = ppe2_m 220 e1v(:,:) = ppe1_m ; e2v(:,:) = ppe2_m 221 e1f(:,:) = ppe1_m ; e2f(:,:) = ppe2_m 222 ! 223 CASE ( 4 ) !== geographical mesh on the sphere, isotropic MERCATOR type ==! 224 ! 225 IF(lwp) WRITE(numout,*) 226 IF(lwp) WRITE(numout,*) ' geographical mesh on the sphere, MERCATOR type' 227 IF(lwp) WRITE(numout,*) ' longitudinal/latitudinal spacing given by ppe1_deg' 228 IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' ) 229 ! 230 ! Find index corresponding to the equator, given the grid spacing e1_deg 231 ! and the (approximate) southern latitude ppgphi0. 232 ! This way we ensure that the equator is at a "T / U" point, when in the domain. 233 ! The formula should work even if the equator is outside the domain. 234 zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2. 235 ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg ) 236 IF( ppgphi0 > 0 ) ijeq = -ijeq 237 ! 238 IF(lwp) WRITE(numout,*) ' Index of the equator on the MERCATOR grid:', ijeq 239 ! 240 DO jj = 1, jpj 241 DO ji = 1, jpi 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 246 ! Longitude 247 glamt(ji,jj) = ppglam0 + ppe1_deg * zti 248 glamu(ji,jj) = ppglam0 + ppe1_deg * zui 249 glamv(ji,jj) = ppglam0 + ppe1_deg * zvi 250 glamf(ji,jj) = ppglam0 + ppe1_deg * zfi 251 ! Latitude 252 gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) ) 253 gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) ) 254 gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) ) 255 gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) ) 256 ! e1 257 e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg 258 e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg 259 e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg 260 e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg 261 ! e2 262 e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg 263 e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg 264 e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg 265 e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg 266 END DO 267 END DO 268 ! 269 CASE ( 5 ) !== beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration) 270 ! 271 IF(lwp) WRITE(numout,*) 272 IF(lwp) WRITE(numout,*) ' beta-plane with regular grid-spacing and rotated domain (GYRE configuration)' 273 IF(lwp) WRITE(numout,*) ' given by ppe1_m and ppe2_m' 274 ! 275 ! Position coordinates (in kilometers) 276 ! ========== 277 ! 278 ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN 279 zlam1 = -85._wp 280 zphi1 = 29._wp 281 ! resolution in meters 282 ze1 = 106000. / REAL( jp_cfg , wp ) 283 ! benchmark: forced the resolution to be about 100 km 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 287 ze1deg = ze1 / (ra * rad) 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 ! 293 IF( nprint==1 .AND. lwp ) THEN 294 WRITE(numout,*) ' ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 295 WRITE(numout,*) ' ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0 296 ENDIF 297 ! 298 DO jj = 1, jpj 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 ! 317 ! Horizontal scale factors (in meters) 318 ! ====== 319 e1t(:,:) = ze1 ; e2t(:,:) = ze1 320 e1u(:,:) = ze1 ; e2u(:,:) = ze1 321 e1v(:,:) = ze1 ; e2v(:,:) = ze1 322 e1f(:,:) = ze1 ; e2f(:,:) = ze1 323 ! 324 CASE DEFAULT 325 WRITE(ctmp1,*) ' bad flag value for jphgr_msh = ', jphgr_msh 326 CALL ctl_stop( ctmp1 ) 327 ! 328 END SELECT 329 330 ! associated horizontal metrics 331 ! ----------------------------- 127 ENDIF 128 129 ! 130 ! !== associated horizontal metrics ==! 332 131 ! 333 132 r1_e1t(:,:) = 1._wp / e1t(:,:) ; r1_e2t (:,:) = 1._wp / e2t(:,:) … … 338 137 e1e2t (:,:) = e1t(:,:) * e2t(:,:) ; r1_e1e2t(:,:) = 1._wp / e1e2t(:,:) 339 138 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(:,:) 139 IF( ie1e2u_v == 0 ) THEN ! u- & v-surfaces have not been defined 140 IF(lwp) WRITE(numout,*) ' u- & v-surfaces calculated as e1 e2 product' 141 e1e2u (:,:) = e1u(:,:) * e2u(:,:) ! compute them 342 142 e1e2v (:,:) = e1v(:,:) * e2v(:,:) 343 ENDIF 344 r1_e1e2u(:,:) = 1._wp / e1e2u(:,:) ! compute their invert in both cases 143 ELSE 144 IF(lwp) WRITE(numout,*) ' u- & v-surfaces have been read in "mesh_mask" file:' 145 IF(lwp) WRITE(numout,*) ' grid size reduction in strait(s) is used' 146 ENDIF 147 r1_e1e2u(:,:) = 1._wp / e1e2u(:,:) ! compute their invert in any cases 345 148 r1_e1e2v(:,:) = 1._wp / e1e2v(:,:) 346 149 ! 347 150 e2_e1u(:,:) = e2u(:,:) / e1u(:,:) 348 151 e1_e2v(:,:) = e1v(:,:) / e2v(:,:) 349 350 IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN ! Control print : Grid informations (if not restart) 351 WRITE(numout,*) 352 WRITE(numout,*) ' longitude and e1 scale factors' 353 WRITE(numout,*) ' ------------------------------' 354 WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1), & 355 glamv(ji,1), glamf(ji,1), & 356 e1t(ji,1), e1u(ji,1), & 357 e1v(ji,1), e1f(ji,1), ji = 1, jpi,10) 358 9300 FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x, & 359 f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 ) 360 ! 361 WRITE(numout,*) 362 WRITE(numout,*) ' latitude and e2 scale factors' 363 WRITE(numout,*) ' -----------------------------' 364 WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj), & 365 & gphiv(1,jj), gphif(1,jj), & 366 & e2t (1,jj), e2u (1,jj), & 367 & e2v (1,jj), e2f (1,jj), jj = 1, jpj, 10 ) 368 ENDIF 369 370 371 ! ================= ! 372 ! Coriolis factor ! 373 ! ================= ! 374 375 SELECT CASE( jphgr_msh ) ! type of horizontal mesh 376 ! 377 CASE ( 0, 1, 4 ) ! mesh on the sphere 378 ! 379 ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 380 ! 381 CASE ( 2 ) ! f-plane at ppgphi0 382 ! 383 ff(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 384 ! 385 IF(lwp) WRITE(numout,*) ' f-plane: Coriolis parameter = constant = ', ff(1,1) 386 ! 387 CASE ( 3 ) ! beta-plane 388 ! 389 zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra ! beta at latitude ppgphi0 390 zphi0 = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad ) ! latitude of the first row F-points 391 ! 392 #if defined key_agrif 393 IF( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN ! for EEL6 configuration only 394 IF( .NOT.Agrif_Root() ) THEN 395 zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad) 396 ENDIF 397 ENDIF 398 #endif 399 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 400 ! 401 ff(:,:) = ( zf0 + zbeta * gphif(:,:) * 1.e+3 ) ! f = f0 +beta* y ( y=0 at south) 402 ! 403 IF(lwp) THEN 404 WRITE(numout,*) 405 WRITE(numout,*) ' Beta-plane: Beta parameter = constant = ', ff(nldi,nldj) 406 WRITE(numout,*) ' Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 407 ENDIF 408 IF( lk_mpp ) THEN 409 zminff=ff(nldi,nldj) 410 zmaxff=ff(nldi,nlej) 411 CALL mpp_min( zminff ) ! min over the global domain 412 CALL mpp_max( zmaxff ) ! max over the global domain 413 IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff 414 END IF 415 ! 416 CASE ( 5 ) ! beta-plane and rotated domain (gyre configuration) 417 ! 418 zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra ! beta at latitude ppgphi0 419 zphi0 = 15._wp ! latitude of the first row F-points 420 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south 421 ! 422 ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south) 423 ! 424 IF(lwp) THEN 425 WRITE(numout,*) 426 WRITE(numout,*) ' Beta-plane and rotated domain : ' 427 WRITE(numout,*) ' Coriolis parameter varies in this processor from ', ff(nldi,nldj),' to ', ff(nldi,nlej) 428 ENDIF 429 ! 430 IF( lk_mpp ) THEN 431 zminff=ff(nldi,nldj) 432 zmaxff=ff(nldi,nlej) 433 CALL mpp_min( zminff ) ! min over the global domain 434 CALL mpp_max( zmaxff ) ! max over the global domain 435 IF(lwp) WRITE(numout,*) ' Coriolis parameter varies globally from ', zminff,' to ', zmaxff 436 END IF 437 ! 438 END SELECT 439 440 441 ! Control of domain for symetrical condition 442 ! ------------------------------------------ 443 ! The equator line must be the latitude coordinate axe 444 445 IF( nperio == 2 ) THEN 446 znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 447 IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' ) 448 ENDIF 152 ! 449 153 ! 450 154 IF( nn_timing == 1 ) CALL timing_stop('dom_hgr') … … 453 157 454 158 455 SUBROUTINE hgr_read( ke1e2u_v ) 159 SUBROUTINE hgr_read( plamt , plamu , plamv , plamf , & ! gridpoints position (required) 160 & pphit , pphiu , pphiv , pphif , & 161 & kff , pff_f , pff_t , & ! Coriolis parameter (if not on the sphere) 162 & pe1t , pe1u , pe1v , pe1f , & ! scale factors (required) 163 & pe2t , pe2u , pe2v , pe2f , & 164 & ke1e2u_v , pe1e2u , pe1e2v ) ! u- & v-surfaces (if gridsize reduction in some straits) 456 165 !!--------------------------------------------------------------------- 457 166 !! *** ROUTINE hgr_read *** 458 167 !! 459 !! ** Purpose : Read a coordinate file in NetCDF format using IOM 460 !! 461 !!---------------------------------------------------------------------- 462 USE iom 463 !! 464 INTEGER, INTENT( inout ) :: ke1e2u_v ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0) 465 ! 466 INTEGER :: inum ! temporary logical unit 168 !! ** Purpose : Read a mesh_mask file in NetCDF format using IOM 169 !! 170 !!---------------------------------------------------------------------- 171 REAL(wp), DIMENSION(:,:), INTENT(out) :: plamt, plamu, plamv, plamf ! longitude outputs 172 REAL(wp), DIMENSION(:,:), INTENT(out) :: pphit, pphiu, pphiv, pphif ! latitude outputs 173 INTEGER , INTENT(out) :: kff ! =1 Coriolis parameter read here, =0 otherwise 174 REAL(wp), DIMENSION(:,:), INTENT(out) :: pff_f, pff_t ! Coriolis factor at f-point (if found in file) 175 REAL(wp), DIMENSION(:,:), INTENT(out) :: pe1t, pe1u, pe1v, pe1f ! i-scale factors 176 REAL(wp), DIMENSION(:,:), INTENT(out) :: pe2t, pe2u, pe2v, pe2f ! j-scale factors 177 INTEGER , INTENT(out) :: ke1e2u_v ! =1 u- & v-surfaces read here, =0 otherwise 178 REAL(wp), DIMENSION(:,:), INTENT(out) :: pe1e2u, pe1e2v ! u- & v-surfaces (if found in file) 179 ! 180 INTEGER :: inum ! logical unit 467 181 !!---------------------------------------------------------------------- 468 182 ! 469 183 IF(lwp) THEN 470 184 WRITE(numout,*) 471 WRITE(numout,*) 'hgr_read : read the horizontal coordinates '185 WRITE(numout,*) 'hgr_read : read the horizontal coordinates in mesh_mask' 472 186 WRITE(numout,*) '~~~~~~~~ jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk 473 187 ENDIF 474 188 ! 475 CALL iom_open( 'coordinates', inum ) 476 ! 477 CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 478 CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr ) 479 CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr ) 480 CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr ) 481 ! 482 CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr ) 483 CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr ) 484 CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr ) 485 CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr ) 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 ) 189 CALL iom_open( cn_domcfg, inum ) 190 ! 191 CALL iom_get( inum, jpdom_data, 'glamt', plamt, lrowattr=ln_use_jattr ) 192 CALL iom_get( inum, jpdom_data, 'glamu', plamu, lrowattr=ln_use_jattr ) 193 CALL iom_get( inum, jpdom_data, 'glamv', plamv, lrowattr=ln_use_jattr ) 194 CALL iom_get( inum, jpdom_data, 'glamf', plamf, lrowattr=ln_use_jattr ) 195 ! 196 CALL iom_get( inum, jpdom_data, 'gphit', pphit, lrowattr=ln_use_jattr ) 197 CALL iom_get( inum, jpdom_data, 'gphiu', pphiu, lrowattr=ln_use_jattr ) 198 CALL iom_get( inum, jpdom_data, 'gphiv', pphiv, lrowattr=ln_use_jattr ) 199 CALL iom_get( inum, jpdom_data, 'gphif', pphif, lrowattr=ln_use_jattr ) 200 ! 201 CALL iom_get( inum, jpdom_data, 'e1t' , pe1t , lrowattr=ln_use_jattr ) 202 CALL iom_get( inum, jpdom_data, 'e1u' , pe1u , lrowattr=ln_use_jattr ) 203 CALL iom_get( inum, jpdom_data, 'e1v' , pe1v , lrowattr=ln_use_jattr ) 204 CALL iom_get( inum, jpdom_data, 'e1f' , pe1f , lrowattr=ln_use_jattr ) 205 ! 206 CALL iom_get( inum, jpdom_data, 'e2t' , pe2t , lrowattr=ln_use_jattr ) 207 CALL iom_get( inum, jpdom_data, 'e2u' , pe2u , lrowattr=ln_use_jattr ) 208 CALL iom_get( inum, jpdom_data, 'e2v' , pe2v , lrowattr=ln_use_jattr ) 209 CALL iom_get( inum, jpdom_data, 'e2f' , pe2f , lrowattr=ln_use_jattr ) 210 ! 211 IF( iom_varid( inum, 'ff_f', ldstop = .FALSE. ) > 0 .AND. & 212 & iom_varid( inum, 'ff_t', ldstop = .FALSE. ) > 0 ) THEN 213 IF(lwp) WRITE(numout,*) ' Coriolis factor at f- and t-points read in ', TRIM( cn_domcfg ), ' file' 214 CALL iom_get( inum, jpdom_data, 'ff_f' , pff_f , lrowattr=ln_use_jattr ) 215 CALL iom_get( inum, jpdom_data, 'ff_t' , pff_t , lrowattr=ln_use_jattr ) 216 kff = 1 217 ELSE 218 kff = 0 219 ENDIF 496 220 ! 497 221 IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN 498 IF(lwp) WRITE(numout,*) ' hgr_read : e1e2u & e1e2v read in coordinatesfile'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 )222 IF(lwp) WRITE(numout,*) ' e1e2u & e1e2v read in ', TRIM( cn_domcfg ), ' file' 223 CALL iom_get( inum, jpdom_data, 'e1e2u' , pe1e2u , lrowattr=ln_use_jattr ) 224 CALL iom_get( inum, jpdom_data, 'e1e2v' , pe1e2v , lrowattr=ln_use_jattr ) 501 225 ke1e2u_v = 1 502 226 ELSE … … 505 229 ! 506 230 CALL iom_close( inum ) 507 508 231 ! 232 END SUBROUTINE hgr_read 509 233 510 234 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
