Changeset 5883 for branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
- Timestamp:
- 2015-11-13T08:01:08+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
r5836 r5883 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(:,:) … … 145 114 CASE DEFAULT ; CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' ) 146 115 END SELECT 147 116 ! 148 117 END SUBROUTINE rot_rep 149 118 … … 155 124 !! ** Purpose : Compute angles between model grid lines and the North direction 156 125 !! 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 126 !! ** Method : sinus and cosinus of the angle between the north-south axe 127 !! and the j-direction at t, u, v and f-points 128 !! dot and cross products are used to obtain cos and sin, resp. 129 !! 130 !! ** Action : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf 131 !!---------------------------------------------------------------------- 132 INTEGER :: ji, jj ! dummy loop indices 133 INTEGER :: ierr ! local integer 134 REAL(wp) :: zlam, zphi ! local scalars 135 REAL(wp) :: zlan, zphh ! - - 136 REAL(wp) :: zxnpt, zynpt, znnpt ! x,y components and norm of the vector: T point to North Pole 137 REAL(wp) :: zxnpu, zynpu, znnpu ! x,y components and norm of the vector: U point to North Pole 138 REAL(wp) :: zxnpv, zynpv, znnpv ! x,y components and norm of the vector: V point to North Pole 139 REAL(wp) :: zxnpf, zynpf, znnpf ! x,y components and norm of the vector: F point to North Pole 140 REAL(wp) :: zxvvt, zyvvt, znvvt ! x,y components and norm of the vector: between V points below and above a T point 141 REAL(wp) :: zxffu, zyffu, znffu ! x,y components and norm of the vector: between F points below and above a U point 142 REAL(wp) :: zxffv, zyffv, znffv ! x,y components and norm of the vector: between F points left and right a V point 143 REAL(wp) :: zxuuf, zyuuf, znuuf ! x,y components and norm of the vector: between U points below and above a F point 144 !!---------------------------------------------------------------------- 145 ! 184 146 ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj), & 185 147 & gsinu(jpi,jpj), gcosu(jpi,jpj), & … … 187 149 & gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr ) 188 150 IF(lk_mpp) CALL mpp_sum( ierr ) 189 IF( ierr /= 0 ) CALL ctl_stop( 'angle: unable to allocate arrays' )190 151 IF( ierr /= 0 ) CALL ctl_stop( 'angle: unable to allocate arrays' ) 152 ! 191 153 ! ============================= ! 192 154 ! Compute the cosinus and sinus ! 193 155 ! ============================= ! 194 156 ! (computation done on the north stereographic polar plane) 195 157 ! 196 158 DO jj = 2, jpjm1 197 159 DO ji = fs_2, jpi ! vector opt. 198 199 ! north pole direction & modulous (at t-point) 200 zlam = glamt(ji,jj) 160 ! 161 zlam = glamt(ji,jj) ! north pole direction & modulous (at t-point) 201 162 zphi = gphit(ji,jj) 202 163 zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 203 164 zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 204 165 znnpt = zxnpt*zxnpt + zynpt*zynpt 205 206 ! north pole direction & modulous (at u-point) 207 zlam = glamu(ji,jj) 166 ! 167 zlam = glamu(ji,jj) ! north pole direction & modulous (at u-point) 208 168 zphi = gphiu(ji,jj) 209 169 zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 210 170 zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 211 171 znnpu = zxnpu*zxnpu + zynpu*zynpu 212 213 ! north pole direction & modulous (at v-point) 214 zlam = glamv(ji,jj) 172 ! 173 zlam = glamv(ji,jj) ! north pole direction & modulous (at v-point) 215 174 zphi = gphiv(ji,jj) 216 175 zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 217 176 zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 218 177 znnpv = zxnpv*zxnpv + zynpv*zynpv 219 220 ! north pole direction & modulous (at f-point) 221 zlam = glamf(ji,jj) 178 ! 179 zlam = glamf(ji,jj) ! north pole direction & modulous (at f-point) 222 180 zphi = gphif(ji,jj) 223 181 zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 224 182 zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 225 183 znnpf = zxnpf*zxnpf + zynpf*zynpf 226 227 ! j-direction: v-point segment direction (around t-point) 228 zlam = glamv(ji,jj ) 184 ! 185 zlam = glamv(ji,jj ) ! j-direction: v-point segment direction (around t-point) 229 186 zphi = gphiv(ji,jj ) 230 187 zlan = glamv(ji,jj-1) … … 236 193 znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt ) ) 237 194 znvvt = MAX( znvvt, 1.e-14 ) 238 239 ! j-direction: f-point segment direction (around u-point) 240 zlam = glamf(ji,jj ) 195 ! 196 zlam = glamf(ji,jj ) ! j-direction: f-point segment direction (around u-point) 241 197 zphi = gphif(ji,jj ) 242 198 zlan = glamf(ji,jj-1) … … 248 204 znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu ) ) 249 205 znffu = MAX( znffu, 1.e-14 ) 250 251 ! i-direction: f-point segment direction (around v-point) 252 zlam = glamf(ji ,jj) 206 ! 207 zlam = glamf(ji ,jj) ! i-direction: f-point segment direction (around v-point) 253 208 zphi = gphif(ji ,jj) 254 209 zlan = glamf(ji-1,jj) … … 260 215 znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv ) ) 261 216 znffv = MAX( znffv, 1.e-14 ) 262 263 ! j-direction: u-point segment direction (around f-point) 264 zlam = glamu(ji,jj+1) 217 ! 218 zlam = glamu(ji,jj+1) ! j-direction: u-point segment direction (around f-point) 265 219 zphi = gphiu(ji,jj+1) 266 220 zlan = glamu(ji,jj ) … … 272 226 znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf ) ) 273 227 znuuf = MAX( znuuf, 1.e-14 ) 274 275 ! cosinus and sinus using scalar and vectorialproducts228 ! 229 ! ! cosinus and sinus using dot and cross products 276 230 gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt 277 231 gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt 278 232 ! 279 233 gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu 280 234 gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu 281 235 ! 282 236 gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf 283 237 gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf 284 285 ! (caution, rotation of 90 degres) 238 ! 286 239 gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv 287 gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv 288 240 gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv ! (caution, rotation of 90 degres) 241 ! 289 242 END DO 290 243 END DO … … 318 271 ! Lateral boundary conditions ! 319 272 ! =========================== ! 320 321 ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 273 ! ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 322 274 CALL lbc_lnk( gcost, 'T', -1. ) ; CALL lbc_lnk( gsint, 'T', -1. ) 323 275 CALL lbc_lnk( gcosu, 'U', -1. ) ; CALL lbc_lnk( gsinu, 'U', -1. ) 324 276 CALL lbc_lnk( gcosv, 'V', -1. ) ; CALL lbc_lnk( gsinv, 'V', -1. ) 325 277 CALL lbc_lnk( gcosf, 'F', -1. ) ; CALL lbc_lnk( gsinf, 'F', -1. ) 326 278 ! 327 279 END SUBROUTINE angle 328 280 329 281 330 SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, & 331 pte, ptn ) 282 SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, pte, ptn ) 332 283 !!---------------------------------------------------------------------- 333 284 !! *** ROUTINE geo2oce *** … … 335 286 !! ** Purpose : 336 287 !! 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 288 !! ** Method : Change a vector from geocentric to east/north 289 !! 347 290 !!---------------------------------------------------------------------- 348 291 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxx, pyy, pzz 349 292 CHARACTER(len=1) , INTENT(in ) :: cgrid 350 293 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pte, ptn 351 ! !294 ! 352 295 REAL(wp), PARAMETER :: rpi = 3.141592653e0 353 296 REAL(wp), PARAMETER :: rad = rpi / 180.e0 … … 355 298 INTEGER :: ierr ! local integer 356 299 !!---------------------------------------------------------------------- 357 300 ! 358 301 IF( .NOT. ALLOCATED( gsinlon ) ) THEN 359 302 ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , & … … 361 304 IF( lk_mpp ) CALL mpp_sum( ierr ) 362 305 IF( ierr /= 0 ) CALL ctl_stop('geo2oce: unable to allocate arrays' ) 306 ENDIF 307 ! 308 SELECT CASE( cgrid) 309 CASE ( 'T' ) 310 ig = 1 311 IF( .NOT. linit(ig) ) THEN 312 gsinlon(:,:,ig) = SIN( rad * glamt(:,:) ) 313 gcoslon(:,:,ig) = COS( rad * glamt(:,:) ) 314 gsinlat(:,:,ig) = SIN( rad * gphit(:,:) ) 315 gcoslat(:,:,ig) = COS( rad * gphit(:,:) ) 316 linit(ig) = .TRUE. 317 ENDIF 318 CASE ( 'U' ) 319 ig = 2 320 IF( .NOT. linit(ig) ) THEN 321 gsinlon(:,:,ig) = SIN( rad * glamu(:,:) ) 322 gcoslon(:,:,ig) = COS( rad * glamu(:,:) ) 323 gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) ) 324 gcoslat(:,:,ig) = COS( rad * gphiu(:,:) ) 325 linit(ig) = .TRUE. 326 ENDIF 327 CASE ( 'V' ) 328 ig = 3 329 IF( .NOT. linit(ig) ) THEN 330 gsinlon(:,:,ig) = SIN( rad * glamv(:,:) ) 331 gcoslon(:,:,ig) = COS( rad * glamv(:,:) ) 332 gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) ) 333 gcoslat(:,:,ig) = COS( rad * gphiv(:,:) ) 334 linit(ig) = .TRUE. 335 ENDIF 336 CASE ( 'F' ) 337 ig = 4 338 IF( .NOT. linit(ig) ) THEN 339 gsinlon(:,:,ig) = SIN( rad * glamf(:,:) ) 340 gcoslon(:,:,ig) = COS( rad * glamf(:,:) ) 341 gsinlat(:,:,ig) = SIN( rad * gphif(:,:) ) 342 gcoslat(:,:,ig) = COS( rad * gphif(:,:) ) 343 linit(ig) = .TRUE. 344 ENDIF 345 CASE default 346 WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid 347 CALL ctl_stop( ctmp1 ) 348 END SELECT 349 ! 350 pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy 351 ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx & 352 & - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy & 353 & + gcoslat(:,:,ig) * pzz 354 ! 355 END SUBROUTINE geo2oce 356 357 358 SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz ) 359 !!---------------------------------------------------------------------- 360 !! *** ROUTINE oce2geo *** 361 !! 362 !! ** Purpose : 363 !! 364 !! ** Method : Change vector from east/north to geocentric 365 !! 366 !! History : ! (A. Caubel) oce2geo - Original code 367 !!---------------------------------------------------------------------- 368 REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pte, ptn 369 CHARACTER(len=1) , INTENT( IN ) :: cgrid 370 REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ) :: pxx , pyy , pzz 371 !! 372 REAL(wp), PARAMETER :: rpi = 3.141592653E0 373 REAL(wp), PARAMETER :: rad = rpi / 180.e0 374 INTEGER :: ig ! 375 INTEGER :: ierr ! local integer 376 !!---------------------------------------------------------------------- 377 378 IF( .NOT. ALLOCATED( gsinlon ) ) THEN 379 ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , & 380 & gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 381 IF( lk_mpp ) CALL mpp_sum( ierr ) 382 IF( ierr /= 0 ) CALL ctl_stop('oce2geo: unable to allocate arrays' ) 363 383 ENDIF 364 384 … … 404 424 CALL ctl_stop( ctmp1 ) 405 425 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 426 ! 427 pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn 428 pyy = gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn 429 pzz = gcoslat(:,:,ig) * ptn 430 ! 493 431 END SUBROUTINE oce2geo 494 432 495 433 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 ) 434 SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv ) 541 435 !!---------------------------------------------------------------------- 542 436 !! *** ROUTINE obs_rot *** … … 546 440 !! current at observation points 547 441 !! 548 !! History : 549 !! 9.2 ! 09-02 (K. Mogensen) 442 !! History : 9.2 ! 09-02 (K. Mogensen) 550 443 !!---------------------------------------------------------------------- 551 444 REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ):: psinu, pcosu, psinv, pcosv ! copy of data 552 445 !!---------------------------------------------------------------------- 553 446 ! 554 447 ! Initialization of gsin* and gcos* at first call 555 448 ! ----------------------------------------------- 556 557 449 IF( lmust_init ) THEN 558 450 IF(lwp) WRITE(numout,*) 559 451 IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched' 560 452 IF(lwp) WRITE(numout,*) ' ~~~~~~~ coordinate transformation' 561 562 453 CALL angle ! initialization of the transformation 563 454 lmust_init = .FALSE. 564 565 455 ENDIF 566 456 ! 567 457 psinu(:,:) = gsinu(:,:) 568 458 pcosu(:,:) = gcosu(:,:) 569 459 psinv(:,:) = gsinv(:,:) 570 460 pcosv(:,:) = gcosv(:,:) 571 461 ! 572 462 END SUBROUTINE obs_rot 573 463
Note: See TracChangeset
for help on using the changeset viewer.