- Timestamp:
- 2016-11-28T17:04:10+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
r5836 r7351 5 5 !!====================================================================== 6 6 !! History : OPA ! 07-1996 (O. Marti) Original code 7 !! NEMO 1.0 ! 02-2008 (G. Madec) F90: Free form 8 !! 3.0 ! 7 !! NEMO 1.0 ! 06-2006 (G. Madec ) Free form, F90 + opt. 8 !! ! 04-2007 (S. Masson) angle: Add T, F points and bugfix in cos lateral boundary 9 !! 3.0 ! 07-2008 (G. Madec) geo2oce suppress lon/lat agruments 10 !! 3.7 ! 11-2015 (G. Madec) remove the unused repere and repcmo routines 9 11 !!---------------------------------------------------------------------- 10 12 11 13 !!---------------------------------------------------------------------- 12 !! r epcmo :13 !! angle :14 !! geo2oce :15 !! repere : old routine suppress it ???14 !! rot_rep : Rotate the Repere: geographic grid <==> stretched coordinates grid 15 !! angle : 16 !! geo2oce : 17 !! oce2geo : 16 18 !!---------------------------------------------------------------------- 17 USE dom_oce ! mesh and scale factors 18 USE phycst ! physical constants 19 USE in_out_manager ! I/O manager 20 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 21 USE lib_mpp ! MPP library 19 USE dom_oce ! mesh and scale factors 20 USE phycst ! physical constants 21 ! 22 USE in_out_manager ! I/O manager 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 24 USE lib_mpp ! MPP library 22 25 23 26 IMPLICIT NONE 24 27 PRIVATE 25 28 26 PUBLIC rot_rep, repcmo, repere, geo2oce, oce2geo ! only rot_rep should be used 27 ! repcmo and repere are keep only for compatibility. 28 ! they are only a useless overlay of rot_rep 29 30 PUBLIC obs_rot 31 32 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: & 33 gsint, gcost, & ! cos/sin between model grid lines and NP direction at T point 34 gsinu, gcosu, & ! cos/sin between model grid lines and NP direction at U point 35 gsinv, gcosv, & ! cos/sin between model grid lines and NP direction at V point 36 gsinf, gcosf ! cos/sin between model grid lines and NP direction at F point 29 PUBLIC rot_rep ! called in sbccpl, fldread, and cyclone 30 PUBLIC geo2oce ! called in sbccpl 31 PUBLIC oce2geo ! called in sbccpl 32 PUBLIC obs_rot ! called in obs_rot_vel and obs_write 33 34 ! ! cos/sin between model grid lines and NP direction 35 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsint, gcost ! at T point 36 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsinu, gcosu ! at U point 37 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsinv, gcosv ! at V point 38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsinf, gcosf ! at F point 37 39 38 40 LOGICAL , SAVE, DIMENSION(4) :: linit = .FALSE. 39 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gsinlon, gcoslon, gsinlat, gcoslat 40 42 41 LOGICAL :: lmust_init = .TRUE. !: used to initialize the cos/sin variables (se above)43 LOGICAL :: lmust_init = .TRUE. !: used to initialize the cos/sin variables (see above) 42 44 43 45 !! * Substitutions … … 50 52 CONTAINS 51 53 52 SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1, &53 px2 , py2 )54 !!----------------------------------------------------------------------55 !! *** ROUTINE repcmo ***56 !!57 !! ** Purpose : Change vector componantes from a geographic grid to a58 !! stretched coordinates grid.59 !!60 !! ** Method : Initialization of arrays at the first call.61 !!62 !! ** Action : - px2 : first componante (defined at u point)63 !! - py2 : second componante (defined at v point)64 !!----------------------------------------------------------------------65 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pxu1, pyu1 ! geographic vector componantes at u-point66 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pxv1, pyv1 ! geographic vector componantes at v-point67 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: px2 ! i-componante (defined at u-point)68 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: py2 ! j-componante (defined at v-point)69 !!----------------------------------------------------------------------70 71 ! Change from geographic to stretched coordinate72 ! ----------------------------------------------73 CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )74 CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )75 76 END SUBROUTINE repcmo77 78 79 54 SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot ) 80 55 !!---------------------------------------------------------------------- … … 83 58 !! ** Purpose : Rotate the Repere: Change vector componantes between 84 59 !! geographic grid <--> stretched coordinates grid. 85 !! 86 !! History : 87 !! 9.2 ! 07-04 (S. Masson) 88 !! (O. Marti ) Original code (repere and repcmo) 89 !!---------------------------------------------------------------------- 90 REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pxin, pyin ! vector componantes 91 CHARACTER(len=1), INTENT( IN ) :: cd_type ! define the nature of pt2d array grid-points 92 CHARACTER(len=5), INTENT( IN ) :: cdtodo ! specify the work to do: 93 !! ! 'en->i' east-north componantes to model i componante 94 !! ! 'en->j' east-north componantes to model j componante 95 !! ! 'ij->e' model i-j componantes to east componante 96 !! ! 'ij->n' model i-j componantes to east componante 97 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: prot 98 !!---------------------------------------------------------------------- 99 100 ! Initialization of gsin* and gcos* at first call 101 ! ----------------------------------------------- 102 103 IF( lmust_init ) THEN 60 !!---------------------------------------------------------------------- 61 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxin, pyin ! vector componantes 62 CHARACTER(len=1), INTENT(in ) :: cd_type ! define the nature of pt2d array grid-points 63 CHARACTER(len=5), INTENT(in ) :: cdtodo ! type of transpormation: 64 ! ! 'en->i' = east-north to i-component 65 ! ! 'en->j' = east-north to j-component 66 ! ! 'ij->e' = (i,j) components to east 67 ! ! 'ij->n' = (i,j) components to north 68 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: prot 69 !!---------------------------------------------------------------------- 70 ! 71 IF( lmust_init ) THEN ! at 1st call only: set gsin. & gcos. 104 72 IF(lwp) WRITE(numout,*) 105 IF(lwp) WRITE(numout,*) ' rot_rep : geographic <--> stretched'106 IF(lwp) WRITE(numout,*) ' ~~~~~ coordinate transformation'73 IF(lwp) WRITE(numout,*) ' rot_rep: coordinate transformation : geographic <==> model (i,j)-components' 74 IF(lwp) WRITE(numout,*) ' ~~~~~~~~ ' 107 75 ! 108 76 CALL angle ! initialization of the transformation 109 77 lmust_init = .FALSE. 110 78 ENDIF 111 112 SELECT CASE (cdtodo) 113 CASE ('en->i') ! 'en->i' est-north componantes to model i componante 79 ! 80 SELECT CASE( cdtodo ) ! type of rotation 81 ! 82 CASE( 'en->i' ) ! east-north to i-component 114 83 SELECT CASE (cd_type) 115 84 CASE ('T') ; prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:) … … 119 88 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 120 89 END SELECT 121 CASE ('en->j') ! 'en->j' est-north componantes to model j componante90 CASE ('en->j') ! east-north to j-component 122 91 SELECT CASE (cd_type) 123 92 CASE ('T') ; prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:) … … 127 96 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 128 97 END SELECT 129 CASE ('ij->e') ! 'ij->e' model i-j componantes to est componante98 CASE ('ij->e') ! (i,j)-components to east 130 99 SELECT CASE (cd_type) 131 100 CASE ('T') ; prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:) … … 135 104 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 136 105 END SELECT 137 CASE ('ij->n') ! 'ij->n' model i-j componantes to est componante106 CASE ('ij->n') ! (i,j)-components to north 138 107 SELECT CASE (cd_type) 139 108 CASE ('T') ; prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:) … … 144 113 END SELECT 145 114 CASE DEFAULT ; CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' ) 115 ! 146 116 END SELECT 147 117 ! 148 118 END SUBROUTINE rot_rep 149 119 … … 155 125 !! ** Purpose : Compute angles between model grid lines and the North direction 156 126 !! 157 !! ** Method : 158 !! 159 !! ** Action : Compute (gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf) arrays: 160 !! sinus and cosinus of the angle between the north-south axe and the 161 !! j-direction at t, u, v and f-points 162 !! 163 !! History : 164 !! 7.0 ! 96-07 (O. Marti ) Original code 165 !! 8.0 ! 98-06 (G. Madec ) 166 !! 8.5 ! 98-06 (G. Madec ) Free form, F90 + opt. 167 !! 9.2 ! 07-04 (S. Masson) Add T, F points and bugfix in cos lateral boundary 168 !!---------------------------------------------------------------------- 169 INTEGER :: ji, jj ! dummy loop indices 170 INTEGER :: ierr ! local integer 171 REAL(wp) :: & 172 zlam, zphi, & ! temporary scalars 173 zlan, zphh, & ! " " 174 zxnpt, zynpt, znnpt, & ! x,y components and norm of the vector: T point to North Pole 175 zxnpu, zynpu, znnpu, & ! x,y components and norm of the vector: U point to North Pole 176 zxnpv, zynpv, znnpv, & ! x,y components and norm of the vector: V point to North Pole 177 zxnpf, zynpf, znnpf, & ! x,y components and norm of the vector: F point to North Pole 178 zxvvt, zyvvt, znvvt, & ! x,y components and norm of the vector: between V points below and above a T point 179 zxffu, zyffu, znffu, & ! x,y components and norm of the vector: between F points below and above a U point 180 zxffv, zyffv, znffv, & ! x,y components and norm of the vector: between F points left and right a V point 181 zxuuf, zyuuf, znuuf ! x,y components and norm of the vector: between U points below and above a F point 182 !!---------------------------------------------------------------------- 183 127 !! ** Method : sinus and cosinus of the angle between the north-south axe 128 !! and the j-direction at t, u, v and f-points 129 !! dot and cross products are used to obtain cos and sin, resp. 130 !! 131 !! ** Action : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf 132 !!---------------------------------------------------------------------- 133 INTEGER :: ji, jj ! dummy loop indices 134 INTEGER :: ierr ! local integer 135 REAL(wp) :: zlam, zphi ! local scalars 136 REAL(wp) :: zlan, zphh ! - - 137 REAL(wp) :: zxnpt, zynpt, znnpt ! x,y components and norm of the vector: T point to North Pole 138 REAL(wp) :: zxnpu, zynpu, znnpu ! x,y components and norm of the vector: U point to North Pole 139 REAL(wp) :: zxnpv, zynpv, znnpv ! x,y components and norm of the vector: V point to North Pole 140 REAL(wp) :: zxnpf, zynpf, znnpf ! x,y components and norm of the vector: F point to North Pole 141 REAL(wp) :: zxvvt, zyvvt, znvvt ! x,y components and norm of the vector: between V points below and above a T point 142 REAL(wp) :: zxffu, zyffu, znffu ! x,y components and norm of the vector: between F points below and above a U point 143 REAL(wp) :: zxffv, zyffv, znffv ! x,y components and norm of the vector: between F points left and right a V point 144 REAL(wp) :: zxuuf, zyuuf, znuuf ! x,y components and norm of the vector: between U points below and above a F point 145 !!---------------------------------------------------------------------- 146 ! 184 147 ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj), & 185 148 & gsinu(jpi,jpj), gcosu(jpi,jpj), & … … 187 150 & gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr ) 188 151 IF(lk_mpp) CALL mpp_sum( ierr ) 189 IF( ierr /= 0 ) CALL ctl_stop( 'angle: unable to allocate arrays' )190 152 IF( ierr /= 0 ) CALL ctl_stop( 'angle: unable to allocate arrays' ) 153 ! 191 154 ! ============================= ! 192 155 ! Compute the cosinus and sinus ! 193 156 ! ============================= ! 194 157 ! (computation done on the north stereographic polar plane) 195 158 ! 196 159 DO jj = 2, jpjm1 197 160 DO ji = fs_2, jpi ! vector opt. 198 199 ! north pole direction & modulous (at t-point) 200 zlam = glamt(ji,jj) 161 ! 162 zlam = glamt(ji,jj) ! north pole direction & modulous (at t-point) 201 163 zphi = gphit(ji,jj) 202 164 zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 203 165 zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 204 166 znnpt = zxnpt*zxnpt + zynpt*zynpt 205 206 ! north pole direction & modulous (at u-point) 207 zlam = glamu(ji,jj) 167 ! 168 zlam = glamu(ji,jj) ! north pole direction & modulous (at u-point) 208 169 zphi = gphiu(ji,jj) 209 170 zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 210 171 zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 211 172 znnpu = zxnpu*zxnpu + zynpu*zynpu 212 213 ! north pole direction & modulous (at v-point) 214 zlam = glamv(ji,jj) 173 ! 174 zlam = glamv(ji,jj) ! north pole direction & modulous (at v-point) 215 175 zphi = gphiv(ji,jj) 216 176 zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 217 177 zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 218 178 znnpv = zxnpv*zxnpv + zynpv*zynpv 219 220 ! north pole direction & modulous (at f-point) 221 zlam = glamf(ji,jj) 179 ! 180 zlam = glamf(ji,jj) ! north pole direction & modulous (at f-point) 222 181 zphi = gphif(ji,jj) 223 182 zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 224 183 zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 225 184 znnpf = zxnpf*zxnpf + zynpf*zynpf 226 227 ! j-direction: v-point segment direction (around t-point) 228 zlam = glamv(ji,jj ) 185 ! 186 zlam = glamv(ji,jj ) ! j-direction: v-point segment direction (around t-point) 229 187 zphi = gphiv(ji,jj ) 230 188 zlan = glamv(ji,jj-1) … … 236 194 znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt ) ) 237 195 znvvt = MAX( znvvt, 1.e-14 ) 238 239 ! j-direction: f-point segment direction (around u-point) 240 zlam = glamf(ji,jj ) 196 ! 197 zlam = glamf(ji,jj ) ! j-direction: f-point segment direction (around u-point) 241 198 zphi = gphif(ji,jj ) 242 199 zlan = glamf(ji,jj-1) … … 248 205 znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu ) ) 249 206 znffu = MAX( znffu, 1.e-14 ) 250 251 ! i-direction: f-point segment direction (around v-point) 252 zlam = glamf(ji ,jj) 207 ! 208 zlam = glamf(ji ,jj) ! i-direction: f-point segment direction (around v-point) 253 209 zphi = gphif(ji ,jj) 254 210 zlan = glamf(ji-1,jj) … … 260 216 znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv ) ) 261 217 znffv = MAX( znffv, 1.e-14 ) 262 263 ! j-direction: u-point segment direction (around f-point) 264 zlam = glamu(ji,jj+1) 218 ! 219 zlam = glamu(ji,jj+1) ! j-direction: u-point segment direction (around f-point) 265 220 zphi = gphiu(ji,jj+1) 266 221 zlan = glamu(ji,jj ) … … 272 227 znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf ) ) 273 228 znuuf = MAX( znuuf, 1.e-14 ) 274 275 ! cosinus and sinus using scalar and vectorialproducts229 ! 230 ! ! cosinus and sinus using dot and cross products 276 231 gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt 277 232 gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt 278 233 ! 279 234 gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu 280 235 gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu 281 236 ! 282 237 gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf 283 238 gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf 284 285 ! (caution, rotation of 90 degres) 239 ! 286 240 gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv 287 gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv 288 241 gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv ! (caution, rotation of 90 degres) 242 ! 289 243 END DO 290 244 END DO … … 318 272 ! Lateral boundary conditions ! 319 273 ! =========================== ! 320 321 ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 274 ! ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 322 275 CALL lbc_lnk( gcost, 'T', -1. ) ; CALL lbc_lnk( gsint, 'T', -1. ) 323 276 CALL lbc_lnk( gcosu, 'U', -1. ) ; CALL lbc_lnk( gsinu, 'U', -1. ) 324 277 CALL lbc_lnk( gcosv, 'V', -1. ) ; CALL lbc_lnk( gsinv, 'V', -1. ) 325 278 CALL lbc_lnk( gcosf, 'F', -1. ) ; CALL lbc_lnk( gsinf, 'F', -1. ) 326 279 ! 327 280 END SUBROUTINE angle 328 281 329 282 330 SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, & 331 pte, ptn ) 283 SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, pte, ptn ) 332 284 !!---------------------------------------------------------------------- 333 285 !! *** ROUTINE geo2oce *** … … 335 287 !! ** Purpose : 336 288 !! 337 !! ** Method : Change wind stress from geocentric to east/north 338 !! 339 !! History : 340 !! ! (O. Marti) Original code 341 !! ! 91-03 (G. Madec) 342 !! ! 92-07 (M. Imbard) 343 !! ! 99-11 (M. Imbard) NetCDF format with IOIPSL 344 !! ! 00-08 (D. Ludicone) Reduced section at Bab el Mandeb 345 !! 8.5 ! 02-06 (G. Madec) F90: Free form 346 !! 3.0 ! 07-08 (G. Madec) geo2oce suppress lon/lat agruments 289 !! ** Method : Change a vector from geocentric to east/north 290 !! 347 291 !!---------------------------------------------------------------------- 348 292 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxx, pyy, pzz 349 293 CHARACTER(len=1) , INTENT(in ) :: cgrid 350 294 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pte, ptn 351 ! !295 ! 352 296 REAL(wp), PARAMETER :: rpi = 3.141592653e0 353 297 REAL(wp), PARAMETER :: rad = rpi / 180.e0 … … 355 299 INTEGER :: ierr ! local integer 356 300 !!---------------------------------------------------------------------- 357 301 ! 358 302 IF( .NOT. ALLOCATED( gsinlon ) ) THEN 359 303 ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , & … … 361 305 IF( lk_mpp ) CALL mpp_sum( ierr ) 362 306 IF( ierr /= 0 ) CALL ctl_stop('geo2oce: unable to allocate arrays' ) 307 ENDIF 308 ! 309 SELECT CASE( cgrid) 310 CASE ( 'T' ) 311 ig = 1 312 IF( .NOT. linit(ig) ) THEN 313 gsinlon(:,:,ig) = SIN( rad * glamt(:,:) ) 314 gcoslon(:,:,ig) = COS( rad * glamt(:,:) ) 315 gsinlat(:,:,ig) = SIN( rad * gphit(:,:) ) 316 gcoslat(:,:,ig) = COS( rad * gphit(:,:) ) 317 linit(ig) = .TRUE. 318 ENDIF 319 CASE ( 'U' ) 320 ig = 2 321 IF( .NOT. linit(ig) ) THEN 322 gsinlon(:,:,ig) = SIN( rad * glamu(:,:) ) 323 gcoslon(:,:,ig) = COS( rad * glamu(:,:) ) 324 gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) ) 325 gcoslat(:,:,ig) = COS( rad * gphiu(:,:) ) 326 linit(ig) = .TRUE. 327 ENDIF 328 CASE ( 'V' ) 329 ig = 3 330 IF( .NOT. linit(ig) ) THEN 331 gsinlon(:,:,ig) = SIN( rad * glamv(:,:) ) 332 gcoslon(:,:,ig) = COS( rad * glamv(:,:) ) 333 gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) ) 334 gcoslat(:,:,ig) = COS( rad * gphiv(:,:) ) 335 linit(ig) = .TRUE. 336 ENDIF 337 CASE ( 'F' ) 338 ig = 4 339 IF( .NOT. linit(ig) ) THEN 340 gsinlon(:,:,ig) = SIN( rad * glamf(:,:) ) 341 gcoslon(:,:,ig) = COS( rad * glamf(:,:) ) 342 gsinlat(:,:,ig) = SIN( rad * gphif(:,:) ) 343 gcoslat(:,:,ig) = COS( rad * gphif(:,:) ) 344 linit(ig) = .TRUE. 345 ENDIF 346 CASE default 347 WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid 348 CALL ctl_stop( ctmp1 ) 349 END SELECT 350 ! 351 pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy 352 ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx & 353 & - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy & 354 & + gcoslat(:,:,ig) * pzz 355 ! 356 END SUBROUTINE geo2oce 357 358 359 SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz ) 360 !!---------------------------------------------------------------------- 361 !! *** ROUTINE oce2geo *** 362 !! 363 !! ** Purpose : 364 !! 365 !! ** Method : Change vector from east/north to geocentric 366 !! 367 !! History : ! (A. Caubel) oce2geo - Original code 368 !!---------------------------------------------------------------------- 369 REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pte, ptn 370 CHARACTER(len=1) , INTENT( IN ) :: cgrid 371 REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ) :: pxx , pyy , pzz 372 !! 373 REAL(wp), PARAMETER :: rpi = 3.141592653E0 374 REAL(wp), PARAMETER :: rad = rpi / 180.e0 375 INTEGER :: ig ! 376 INTEGER :: ierr ! local integer 377 !!---------------------------------------------------------------------- 378 379 IF( .NOT. ALLOCATED( gsinlon ) ) THEN 380 ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , & 381 & gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 382 IF( lk_mpp ) CALL mpp_sum( ierr ) 383 IF( ierr /= 0 ) CALL ctl_stop('oce2geo: unable to allocate arrays' ) 363 384 ENDIF 364 385 … … 404 425 CALL ctl_stop( ctmp1 ) 405 426 END SELECT 406 407 pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy 408 ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx & 409 - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy & 410 + gcoslat(:,:,ig) * pzz 411 !!$ ptv = gcoslon(:,:,ig) * gcoslat(:,:,ig) * pxx & 412 !!$ + gsinlon(:,:,ig) * gcoslat(:,:,ig) * pyy & 413 !!$ + gsinlat(:,:,ig) * pzz 414 ! 415 END SUBROUTINE geo2oce 416 417 SUBROUTINE oce2geo ( pte, ptn, cgrid, & 418 pxx , pyy , pzz ) 419 !!---------------------------------------------------------------------- 420 !! *** ROUTINE oce2geo *** 421 !! 422 !! ** Purpose : 423 !! 424 !! ** Method : Change vector from east/north to geocentric 425 !! 426 !! History : 427 !! ! (A. Caubel) oce2geo - Original code 428 !!---------------------------------------------------------------------- 429 REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pte, ptn 430 CHARACTER(len=1) , INTENT( IN ) :: cgrid 431 REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ) :: pxx , pyy , pzz 432 !! 433 REAL(wp), PARAMETER :: rpi = 3.141592653E0 434 REAL(wp), PARAMETER :: rad = rpi / 180.e0 435 INTEGER :: ig ! 436 INTEGER :: ierr ! local integer 437 !!---------------------------------------------------------------------- 438 439 IF( .NOT. ALLOCATED( gsinlon ) ) THEN 440 ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , & 441 & gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 442 IF( lk_mpp ) CALL mpp_sum( ierr ) 443 IF( ierr /= 0 ) CALL ctl_stop('oce2geo: unable to allocate arrays' ) 444 ENDIF 445 446 SELECT CASE( cgrid) 447 CASE ( 'T' ) 448 ig = 1 449 IF( .NOT. linit(ig) ) THEN 450 gsinlon(:,:,ig) = SIN( rad * glamt(:,:) ) 451 gcoslon(:,:,ig) = COS( rad * glamt(:,:) ) 452 gsinlat(:,:,ig) = SIN( rad * gphit(:,:) ) 453 gcoslat(:,:,ig) = COS( rad * gphit(:,:) ) 454 linit(ig) = .TRUE. 455 ENDIF 456 CASE ( 'U' ) 457 ig = 2 458 IF( .NOT. linit(ig) ) THEN 459 gsinlon(:,:,ig) = SIN( rad * glamu(:,:) ) 460 gcoslon(:,:,ig) = COS( rad * glamu(:,:) ) 461 gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) ) 462 gcoslat(:,:,ig) = COS( rad * gphiu(:,:) ) 463 linit(ig) = .TRUE. 464 ENDIF 465 CASE ( 'V' ) 466 ig = 3 467 IF( .NOT. linit(ig) ) THEN 468 gsinlon(:,:,ig) = SIN( rad * glamv(:,:) ) 469 gcoslon(:,:,ig) = COS( rad * glamv(:,:) ) 470 gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) ) 471 gcoslat(:,:,ig) = COS( rad * gphiv(:,:) ) 472 linit(ig) = .TRUE. 473 ENDIF 474 CASE ( 'F' ) 475 ig = 4 476 IF( .NOT. linit(ig) ) THEN 477 gsinlon(:,:,ig) = SIN( rad * glamf(:,:) ) 478 gcoslon(:,:,ig) = COS( rad * glamf(:,:) ) 479 gsinlat(:,:,ig) = SIN( rad * gphif(:,:) ) 480 gcoslat(:,:,ig) = COS( rad * gphif(:,:) ) 481 linit(ig) = .TRUE. 482 ENDIF 483 CASE default 484 WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid 485 CALL ctl_stop( ctmp1 ) 486 END SELECT 487 488 pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn 489 pyy = gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn 490 pzz = gcoslat(:,:,ig) * ptn 491 492 427 ! 428 pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn 429 pyy = gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn 430 pzz = gcoslat(:,:,ig) * ptn 431 ! 493 432 END SUBROUTINE oce2geo 494 433 495 434 496 SUBROUTINE repere ( px1, py1, px2, py2, kchoix, cd_type ) 497 !!---------------------------------------------------------------------- 498 !! *** ROUTINE repere *** 499 !! 500 !! ** Purpose : Change vector componantes between a geopgraphic grid 501 !! and a stretched coordinates grid. 502 !! 503 !! ** Method : 504 !! 505 !! ** Action : 506 !! 507 !! History : 508 !! ! 89-03 (O. Marti) original code 509 !! ! 92-02 (M. Imbard) 510 !! ! 93-03 (M. Guyon) symetrical conditions 511 !! ! 98-05 (B. Blanke) 512 !! 8.5 ! 02-08 (G. Madec) F90: Free form 513 !!---------------------------------------------------------------------- 514 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: px1, py1 ! two horizontal components to be rotated 515 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: px2, py2 ! the two horizontal components in the model repere 516 INTEGER , INTENT(in ) :: kchoix ! type of transformation 517 ! ! = 1 change from geographic to model grid. 518 ! ! =-1 change from model to geographic grid 519 CHARACTER(len=1), INTENT(in ), OPTIONAL :: cd_type ! define the nature of pt2d array grid-points 520 ! 521 CHARACTER(len=1) :: cl_type ! define the nature of pt2d array grid-points (T point by default) 522 !!---------------------------------------------------------------------- 523 524 cl_type = 'T' 525 IF( PRESENT(cd_type) ) cl_type = cd_type 526 ! 527 SELECT CASE (kchoix) 528 CASE ( 1) ! change from geographic to model grid. 529 CALL rot_rep( px1, py1, cl_type, 'en->i', px2 ) 530 CALL rot_rep( px1, py1, cl_type, 'en->j', py2 ) 531 CASE (-1) ! change from model to geographic grid 532 CALL rot_rep( px1, py1, cl_type, 'ij->e', px2 ) 533 CALL rot_rep( px1, py1, cl_type, 'ij->n', py2 ) 534 CASE DEFAULT ; CALL ctl_stop( 'repere: Syntax Error in the definition of kchoix (1 OR -1' ) 535 END SELECT 536 537 END SUBROUTINE repere 538 539 540 SUBROUTINE obs_rot ( psinu, pcosu, psinv, pcosv ) 435 SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv ) 541 436 !!---------------------------------------------------------------------- 542 437 !! *** ROUTINE obs_rot *** … … 546 441 !! current at observation points 547 442 !! 548 !! History : 549 !! 9.2 ! 09-02 (K. Mogensen) 443 !! History : 9.2 ! 09-02 (K. Mogensen) 550 444 !!---------------------------------------------------------------------- 551 445 REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ):: psinu, pcosu, psinv, pcosv ! copy of data 552 446 !!---------------------------------------------------------------------- 553 447 ! 554 448 ! Initialization of gsin* and gcos* at first call 555 449 ! ----------------------------------------------- 556 557 450 IF( lmust_init ) THEN 558 451 IF(lwp) WRITE(numout,*) 559 452 IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched' 560 453 IF(lwp) WRITE(numout,*) ' ~~~~~~~ coordinate transformation' 561 562 454 CALL angle ! initialization of the transformation 563 455 lmust_init = .FALSE. 564 565 456 ENDIF 566 457 ! 567 458 psinu(:,:) = gsinu(:,:) 568 459 pcosu(:,:) = gcosu(:,:) 569 460 psinv(:,:) = gsinv(:,:) 570 461 pcosv(:,:) = gcosv(:,:) 571 462 ! 572 463 END SUBROUTINE obs_rot 573 464
Note: See TracChangeset
for help on using the changeset viewer.