- Timestamp:
- 2011-02-27T17:33:43+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
r2620 r2624 4 4 !! Ocean mesh : ??? 5 5 !!====================================================================== 6 !! History : OPA ! 1989-07 (O. Marti) Original code 7 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form + opt. 8 !! - ! 2006-02 (A. Caubel) oce2geo - Original code 9 !! 2.0 ! 2007-04 (S. Masson) Add T, F points and bugfix in cos lateral boundary 10 !! 3.0 ! 2008-07 (G. Madec) geo2oce suppress lon/lat agruments 11 !! 3.3 ! 2010-10 (K. Mogensen) add obs_rot 12 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation + addition of angle_geo 6 !! History : OPA ! 07-1996 (O. Marti) Original code 7 !! NEMO 1.0 ! 02-2008 (G. Madec) F90: Free form 8 !! 3.0 ! 13 9 !!---------------------------------------------------------------------- 14 10 15 11 !!---------------------------------------------------------------------- 16 !! angle_msh_geo : local sin/cos of the angle between model grid lines and the North direction 17 !! angle_geo : local sin/cos of the latitude and longitude at each mesh grid point 18 !! geo2oce, oce2geo : ? 19 !! obs_rot : provide to OBS operator the sin &cos at u- and v-points 20 !! 21 !! repcmo, repere, rot_rep : old routines ==> to be suppressed 12 !! repcmo : 13 !! angle : 14 !! geo2oce : 15 !! repere : old routine suppress it ??? 22 16 !!---------------------------------------------------------------------- 23 USE dom_oce ! mesh and scale factors 24 USE phycst ! physical constants 25 USE in_out_manager ! I/O manager 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 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 27 22 28 23 IMPLICIT NONE 29 24 PRIVATE 30 25 31 PUBLIC geo2oce, oce2geo ! 32 PUBLIC obs_rot ! 33 ! 34 PUBLIC repcmo, repere, rot_rep ! CAUTION: these routines are kept only for compatibility. 35 ! ! They are only a useless overlay of rot_rep 36 37 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gsin, gcos ! sinus & cosinus at T,U,V,F points 38 ! ! between mesh and NP direction 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gsinlon, gsinlat ! sinus of lon & lat at T,U,V,F points 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gcoslon, gcoslat ! cosinus of lon & lat at T,U,V,F points 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 37 38 LOGICAL , SAVE, DIMENSION(4) :: linit = .FALSE. 39 REAL(wp), SAVE, DIMENSION(:,:,:) :: gsinlon, gcoslon, gsinlat, gcoslat 40 41 LOGICAL :: lmust_init = .TRUE. !: used to initialize the cos/sin variables (se above) 41 42 42 43 !! * Substitutions 43 44 # include "vectopt_loop_substitute.h90" 44 45 !!---------------------------------------------------------------------- 45 !! NEMO/OPA 4.0 , NEMO Consortium (2011)46 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 46 47 !! $Id$ 47 !! Software governed by the CeCILL licence 48 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 48 49 !!---------------------------------------------------------------------- 50 49 51 CONTAINS 50 52 51 SUBROUTINE angle_msh_geo 52 !!---------------------------------------------------------------------- 53 !! *** ROUTINE angle_msh_geo *** 54 !! 55 !! ** Purpose : Compute angles between model grid lines and the North direction 56 !! 57 !! ** Action : allocate and compute (gsint, gcost, gsinu, gcosu, gsinv, 58 !! gcosv, gsinf, gcosf) arrays: sinus and cosinus of the angle 59 !! between the north-south axe and the j-direction at t, u, v and f-points 60 !!---------------------------------------------------------------------- 61 INTEGER :: ji, jj ! dummy loop indices 62 INTEGER :: ierr ! local integer 63 REAL(wp) :: zpi_4 ! local scalar (pi/4) 64 REAL(wp) :: zlam, zphi, zlan, zphh ! local scalars 65 REAL(wp) :: zxnpt, zynpt, znnpt ! x,y components & norm of the vector: T point to North Pole 66 REAL(wp) :: zxnpu, zynpu, znnpu ! x,y components & norm of the vector: U point to North Pole 67 REAL(wp) :: zxnpv, zynpv, znnpv ! x,y components & norm of the vector: V point to North Pole 68 REAL(wp) :: zxnpf, zynpf, znnpf ! x,y components & norm of the vector: F point to North Pole 69 REAL(wp) :: zxvvt, zyvvt, znvvt ! x,y components & norm of the vector: between V pts below and above a T pt 70 REAL(wp) :: zxffu, zyffu, znffu ! x,y components & norm of the vector: between F pts below and above a U pt 71 REAL(wp) :: zxffv, zyffv, znffv ! x,y components & norm of the vector: between F pts left and right a V pt 72 REAL(wp) :: zxuuf, zyuuf, znuuf ! x,y components & norm of the vector: between U pts below and above a F pt 73 !!---------------------------------------------------------------------- 74 ! 75 ! already allocated & initialized ==> return 76 ! ------------------------------- 77 IF( ALLOCATED( gsin ) .AND. ALLOCATED( gcos ) ) RETURN 78 79 ! allocate cos & sin arrays 80 ! ------------------------- 81 ALLOCATE( gsin(jpi,jpj,4) , gcos(jpi,jpj,4) , STAT=ierr ) 82 IF(lk_mpp) CALL mpp_sum( ierr ) 83 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'angle_msh_geo: unable to allocate arrays' ) 84 ! 85 ! initialize cos & sin arrays 86 ! ---------------------------- 87 ! (computation done on the north stereographic polar plane) 88 IF(lwp) WRITE(numout,*) 89 IF(lwp) WRITE(numout,*) ' angle_msh_geo : geographic <--> model mesh , cos/sin initialization ' 90 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 91 92 zpi_4 = rpi / 4._wp 93 ! 94 DO jj = 2, jpjm1 95 !CDIR NOVERRCHK 96 DO ji = fs_2, jpi ! vector opt. 97 ! 98 zlam = glamt(ji,jj) ! north pole direction & modulous (at t-point) 99 zphi = gphit(ji,jj) 100 zxnpt = 0._wp - 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) 101 zynpt = 0._wp - 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) 102 znnpt = zxnpt*zxnpt + zynpt*zynpt 103 104 zlam = glamu(ji,jj) ! north pole direction & modulous (at u-point) 105 zphi = gphiu(ji,jj) 106 zxnpu = 0._wp - 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) 107 zynpu = 0._wp - 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) 108 znnpu = zxnpu*zxnpu + zynpu*zynpu 109 110 zlam = glamv(ji,jj) ! north pole direction & modulous (at v-point) 111 zphi = gphiv(ji,jj) 112 zxnpv = 0._wp - 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) 113 zynpv = 0._wp - 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) 114 znnpv = zxnpv*zxnpv + zynpv*zynpv 115 116 zlam = glamf(ji,jj) ! north pole direction & modulous (at f-point) 117 zphi = gphif(ji,jj) 118 zxnpf = 0._wp - 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) 119 zynpf = 0._wp - 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) 120 znnpf = zxnpf*zxnpf + zynpf*zynpf 121 122 zlam = glamv(ji,jj ) ! j-direction: v-point segment direction (around t-point) 123 zphi = gphiv(ji,jj ) 124 zlan = glamv(ji,jj-1) 125 zphh = gphiv(ji,jj-1) 126 zxvvt = 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) & 127 & - 2._wp * COS( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 ) 128 zyvvt = 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) & 129 & - 2._wp * SIN( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 ) 130 znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt ) ) 131 znvvt = MAX( znvvt, 1.e-14 ) 132 133 zlam = glamf(ji,jj ) ! j-direction: f-point segment direction (around u-point) 134 zphi = gphif(ji,jj ) 135 zlan = glamf(ji,jj-1) 136 zphh = gphif(ji,jj-1) 137 zxffu = 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) & 138 & - 2._wp * COS( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 ) 139 zyffu = 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) & 140 & - 2._wp * SIN( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 ) 141 znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu ) ) 142 znffu = MAX( znffu, 1.e-14 ) 143 144 zlam = glamf(ji ,jj) ! i-direction: f-point segment direction (around v-point) 145 zphi = gphif(ji ,jj) 146 zlan = glamf(ji-1,jj) 147 zphh = gphif(ji-1,jj) 148 zxffv = 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) & 149 & - 2._wp * COS( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 ) 150 zyffv = 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) & 151 & - 2._wp * SIN( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 ) 152 znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv ) ) 153 znffv = MAX( znffv, 1.e-14 ) 154 155 zlam = glamu(ji,jj+1) ! j-direction: u-point segment direction (around f-point) 156 zphi = gphiu(ji,jj+1) 157 zlan = glamu(ji,jj ) 158 zphh = gphiu(ji,jj ) 159 zxuuf = 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) & 160 & - 2._wp * COS( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 ) 161 zyuuf = 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 ) & 162 & - 2._wp * SIN( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 ) 163 znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf ) ) 164 znuuf = MAX( znuuf, 1.e-14 ) 165 166 ! cosinus and sinus using scalar and vectorial products 167 gsin(ji,jj,1) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt ! T-point 168 gcos(ji,jj,1) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt 169 170 gsin(ji,jj,2) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu ! U-point 171 gcos(ji,jj,2) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu 172 ! 173 gsin(ji,jj,3) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv ! F-point 174 gcos(ji,jj,3) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv ! (caution, rotation of 90 degres) 175 ! 176 gsin(ji,jj,4) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf ! V-point 177 gcos(ji,jj,4) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf 178 ! 179 END DO 180 END DO 181 ! 182 ! Geographic mesh case 183 ! -------------------- 184 DO jj = 2, jpjm1 185 DO ji = fs_2, jpi ! vector opt. 186 IF( MOD( ABS( glamv(ji,jj) - glamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN 187 gsin(ji,jj,1) = 0. 188 gcos(ji,jj,1) = 1. 189 ENDIF 190 IF( MOD( ABS( glamf(ji,jj) - glamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN 191 gsin(ji,jj,2) = 0. 192 gcos(ji,jj,2) = 1. 193 ENDIF 194 IF( ABS( gphif(ji,jj) - gphif(ji-1,jj) ) < 1.e-8 ) THEN 195 gsin(ji,jj,3) = 0. 196 gcos(ji,jj,3) = 1. 197 ENDIF 198 IF( MOD( ABS( glamu(ji,jj) - glamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN 199 gsin(ji,jj,4) = 0. 200 gcos(ji,jj,4) = 1. 201 ENDIF 202 END DO 203 END DO 204 ! 205 CALL lbc_lnk( gcos(:,:,1), 'T', -1. ) ; CALL lbc_lnk( gsin(:,:,1), 'T', -1. ) ! lateral boundary cond. 206 CALL lbc_lnk( gcos(:,:,2), 'U', -1. ) ; CALL lbc_lnk( gsin(:,:,2), 'U', -1. ) 207 CALL lbc_lnk( gcos(:,:,3), 'V', -1. ) ; CALL lbc_lnk( gsin(:,:,3), 'V', -1. ) 208 CALL lbc_lnk( gcos(:,:,4), 'F', -1. ) ; CALL lbc_lnk( gsin(:,:,4), 'F', -1. ) 209 ! 210 END SUBROUTINE angle_msh_geo 211 212 213 SUBROUTINE angle_geo 214 !!---------------------------------------------------------------------- 215 !! *** ROUTINE angle_geo *** 216 !! 217 !! ** Purpose : compute one for all, and at each mesh grid points, 218 !! the local cos/sin of latitude/longitude. 219 !!---------------------------------------------------------------------- 220 INTEGER ierr 221 !!---------------------------------------------------------------------- 222 ! 223 IF( ALLOCATED( gsinlon ) .AND. ALLOCATED( gcoslon ) .AND. & !== already allocated & initialized 224 ALLOCATED( gsinlat ) .AND. ALLOCATED( gcoslat ) ) RETURN 225 226 ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , & !== allocate the arrays 227 & gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 228 IF( lk_mpp ) CALL mpp_sum( ierr ) 229 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'angle_geo: unable to allocate arrays') 230 ! 231 gsinlon(:,:,1) = SIN( rad * glamt(:,:) ) ! T-point 232 gcoslon(:,:,1) = COS( rad * glamt(:,:) ) 233 gsinlat(:,:,1) = SIN( rad * gphit(:,:) ) 234 gcoslat(:,:,1) = COS( rad * gphit(:,:) ) 235 ! 236 gsinlon(:,:,2) = SIN( rad * glamu(:,:) ) ! U-point 237 gcoslon(:,:,2) = COS( rad * glamu(:,:) ) 238 gsinlat(:,:,2) = SIN( rad * gphiu(:,:) ) 239 gcoslat(:,:,2) = COS( rad * gphiu(:,:) ) 240 ! 241 gsinlon(:,:,3) = SIN( rad * glamv(:,:) ) ! V-point 242 gcoslon(:,:,3) = COS( rad * glamv(:,:) ) 243 gsinlat(:,:,3) = SIN( rad * gphiv(:,:) ) 244 gcoslat(:,:,3) = COS( rad * gphiv(:,:) ) 245 ! 246 gsinlon(:,:,4) = SIN( rad * glamf(:,:) ) ! T-point 247 gcoslon(:,:,4) = COS( rad * glamf(:,:) ) 248 gsinlat(:,:,4) = SIN( rad * gphif(:,:) ) 249 gcoslat(:,:,4) = COS( rad * gphif(:,:) ) 250 ! 251 END SUBROUTINE angle_geo 252 253 254 SUBROUTINE geo2oce( pxx, pyy, pzz, cgrid, pte, ptn ) 255 !!---------------------------------------------------------------------- 256 !! *** ROUTINE geo2oce *** 257 !! 258 !! ** Purpose : Change a vector from geocentric to east/north 259 !!---------------------------------------------------------------------- 260 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxx, pyy, pzz ! 261 CHARACTER(len=1) , INTENT(in ) :: cgrid ! type of grid point 262 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pte, ptn ! 263 !! 264 INTEGER :: ipt ! local integer 265 !!---------------------------------------------------------------------- 266 267 CALL angle_geo ! initialization of coc & sin (just a return if not the 1st call) 268 ! 269 SELECT CASE (cgrid) ! type of point 270 CASE ('T') ; ipt = 1 271 CASE ('U') ; ipt = 2 272 CASE ('V') ; ipt = 3 273 CASE ('F') ; ipt = 4 274 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 275 END SELECT 276 277 ! ! transformation 278 pte(:,:) = - gsinlon(:,:,ipt) * pxx(:,:) + gcoslon(:,:,ipt) * pyy(:,:) 279 ptn(:,:) = - gsinlat(:,:,ipt) * gsinlat(:,:,ipt) * pxx(:,:) & 280 & - gsinlon(:,:,ipt) * gsinlat(:,:,ipt) * pyy(:,:) + gcoslat(:,:,ipt) * pzz(:,:) 281 !!$ ptv(:,:) = gsinlat(:,:,ipt) * gcoslat(:,:,ipt) * pxx(:,:) & 282 !!$ + zsinlon(:,:,ipt) * gcoslat(:,:,ipt) * pyy(:,:) + gsinlat(:,:,ipt) * pzz(:,:) 283 ! 284 END SUBROUTINE geo2oce 285 286 287 SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz ) 288 !!---------------------------------------------------------------------- 289 !! *** ROUTINE oce2geo *** 290 !! 291 !! ** Purpose : Change vector from east/north to geocentric 292 !!---------------------------------------------------------------------- 293 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pte, ptn ! 294 CHARACTER(len=1) , INTENT(in ) :: cgrid ! 295 REAL(wp), DIMENSION(:,:), INTENT( out) :: pxx , pyy , pzz ! 296 !! 297 INTEGER :: ipt ! 298 !!---------------------------------------------------------------------- 299 ! 300 CALL angle_geo ! initialization of coc & sin (just a return if not the 1st call) 301 ! 302 SELECT CASE (cgrid) ! type of point 303 CASE ('T') ; ipt = 1 304 CASE ('U') ; ipt = 2 305 CASE ('V') ; ipt = 3 306 CASE ('F') ; ipt = 4 307 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 308 END SELECT 309 ! 310 ! ! transformation 311 pxx(:,:) = - gsinlon(:,:,ipt) * pte(:,:) - gcoslon(:,:,ipt) * gsinlat(:,:,ipt) * ptn(:,:) 312 pyy(:,:) = gcoslon(:,:,ipt) * pte(:,:) - gsinlon(:,:,ipt) * gsinlat(:,:,ipt) * ptn(:,:) 313 pzz(:,:) = gcoslat(:,:,ipt) * ptn(:,:) 314 ! 315 END SUBROUTINE oce2geo 316 317 318 SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv ) 319 !!---------------------------------------------------------------------- 320 !! *** ROUTINE obs_rot *** 321 !! 322 !! ** Purpose : provide to OBS operator the sin &cos at u- and v-points 323 !! in order to rotate currents at observation points 324 !!---------------------------------------------------------------------- 325 REAL(wp), DIMENSION(jpi,jpj), INTENT(out):: psinu, pcosu, psinv, pcosv ! copy of data 326 !!---------------------------------------------------------------------- 327 ! 328 CALL angle_msh_geo ! initialization of coc & sin (just a return if not the 1st call) 329 ! 330 psinu(:,:) = gsin(:,:,2) ! U-point 331 pcosu(:,:) = gcos(:,:,2) 332 psinv(:,:) = gsin(:,:,3) ! V-point 333 pcosv(:,:) = gcos(:,:,3) 334 ! 335 END SUBROUTINE obs_rot 336 337 338 SUBROUTINE repcmo( pxu1, pyu1, pxv1, pyv1, px2 , py2 ) 53 SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1, & 54 px2 , py2 ) 339 55 !!---------------------------------------------------------------------- 340 56 !! *** ROUTINE repcmo *** 341 57 !! 342 58 !! ** Purpose : Change vector componantes from a geographic grid to a 343 !! 59 !! stretched coordinates grid. 344 60 !! 345 61 !! ** Method : Initialization of arrays at the first call. … … 353 69 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: py2 ! j-componante (defined at v-point) 354 70 !!---------------------------------------------------------------------- 355 ! 356 CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 ) ! Change from geographic to stretched coordinate 71 72 ! Change from geographic to stretched coordinate 73 ! ---------------------------------------------- 74 CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 ) 357 75 CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 ) 76 77 END SUBROUTINE repcmo 78 79 80 SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot ) 81 !!---------------------------------------------------------------------- 82 !! *** ROUTINE rot_rep *** 83 !! 84 !! ** Purpose : Rotate the Repere: Change vector componantes between 85 !! geographic grid <--> stretched coordinates grid. 86 !! 87 !! History : 88 !! 9.2 ! 07-04 (S. Masson) 89 !! (O. Marti ) Original code (repere and repcmo) 90 !!---------------------------------------------------------------------- 91 REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pxin, pyin ! vector componantes 92 CHARACTER(len=1), INTENT( IN ) :: cd_type ! define the nature of pt2d array grid-points 93 CHARACTER(len=5), INTENT( IN ) :: cdtodo ! specify the work to do: 94 !! ! 'en->i' east-north componantes to model i componante 95 !! ! 'en->j' east-north componantes to model j componante 96 !! ! 'ij->e' model i-j componantes to east componante 97 !! ! 'ij->n' model i-j componantes to east componante 98 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: prot 99 !!---------------------------------------------------------------------- 100 101 ! Initialization of gsin* and gcos* at first call 102 ! ----------------------------------------------- 103 104 IF( lmust_init ) THEN 105 IF(lwp) WRITE(numout,*) 106 IF(lwp) WRITE(numout,*) ' rot_rep : geographic <--> stretched' 107 IF(lwp) WRITE(numout,*) ' ~~~~~ coordinate transformation' 108 ! 109 CALL angle ! initialization of the transformation 110 lmust_init = .FALSE. 111 ENDIF 112 113 SELECT CASE (cdtodo) 114 CASE ('en->i') ! 'en->i' est-north componantes to model i componante 115 SELECT CASE (cd_type) 116 CASE ('T') ; prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:) 117 CASE ('U') ; prot(:,:) = pxin(:,:) * gcosu(:,:) + pyin(:,:) * gsinu(:,:) 118 CASE ('V') ; prot(:,:) = pxin(:,:) * gcosv(:,:) + pyin(:,:) * gsinv(:,:) 119 CASE ('F') ; prot(:,:) = pxin(:,:) * gcosf(:,:) + pyin(:,:) * gsinf(:,:) 120 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 121 END SELECT 122 CASE ('en->j') ! 'en->j' est-north componantes to model j componante 123 SELECT CASE (cd_type) 124 CASE ('T') ; prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:) 125 CASE ('U') ; prot(:,:) = pyin(:,:) * gcosu(:,:) - pxin(:,:) * gsinu(:,:) 126 CASE ('V') ; prot(:,:) = pyin(:,:) * gcosv(:,:) - pxin(:,:) * gsinv(:,:) 127 CASE ('F') ; prot(:,:) = pyin(:,:) * gcosf(:,:) - pxin(:,:) * gsinf(:,:) 128 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 129 END SELECT 130 CASE ('ij->e') ! 'ij->e' model i-j componantes to est componante 131 SELECT CASE (cd_type) 132 CASE ('T') ; prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:) 133 CASE ('U') ; prot(:,:) = pxin(:,:) * gcosu(:,:) - pyin(:,:) * gsinu(:,:) 134 CASE ('V') ; prot(:,:) = pxin(:,:) * gcosv(:,:) - pyin(:,:) * gsinv(:,:) 135 CASE ('F') ; prot(:,:) = pxin(:,:) * gcosf(:,:) - pyin(:,:) * gsinf(:,:) 136 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 137 END SELECT 138 CASE ('ij->n') ! 'ij->n' model i-j componantes to est componante 139 SELECT CASE (cd_type) 140 CASE ('T') ; prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:) 141 CASE ('U') ; prot(:,:) = pyin(:,:) * gcosu(:,:) + pxin(:,:) * gsinu(:,:) 142 CASE ('V') ; prot(:,:) = pyin(:,:) * gcosv(:,:) + pxin(:,:) * gsinv(:,:) 143 CASE ('F') ; prot(:,:) = pyin(:,:) * gcosf(:,:) + pxin(:,:) * gsinf(:,:) 144 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 145 END SELECT 146 CASE DEFAULT ; CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' ) 147 END SELECT 148 149 END SUBROUTINE rot_rep 150 151 152 SUBROUTINE angle 153 !!---------------------------------------------------------------------- 154 !! *** ROUTINE angle *** 155 !! 156 !! ** Purpose : Compute angles between model grid lines and the North direction 157 !! 158 !! ** Method : 159 !! 160 !! ** Action : Compute (gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf) arrays: 161 !! sinus and cosinus of the angle between the north-south axe and the 162 !! j-direction at t, u, v and f-points 163 !! 164 !! History : 165 !! 7.0 ! 96-07 (O. Marti ) Original code 166 !! 8.0 ! 98-06 (G. Madec ) 167 !! 8.5 ! 98-06 (G. Madec ) Free form, F90 + opt. 168 !! 9.2 ! 07-04 (S. Masson) Add T, F points and bugfix in cos lateral boundary 169 !!---------------------------------------------------------------------- 170 INTEGER :: ji, jj ! dummy loop indices 171 INTEGER :: ierr ! local integer 172 REAL(wp) :: & 173 zlam, zphi, & ! temporary scalars 174 zlan, zphh, & ! " " 175 zxnpt, zynpt, znnpt, & ! x,y components and norm of the vector: T point to North Pole 176 zxnpu, zynpu, znnpu, & ! x,y components and norm of the vector: U point to North Pole 177 zxnpv, zynpv, znnpv, & ! x,y components and norm of the vector: V point to North Pole 178 zxnpf, zynpf, znnpf, & ! x,y components and norm of the vector: F point to North Pole 179 zxvvt, zyvvt, znvvt, & ! x,y components and norm of the vector: between V points below and above a T point 180 zxffu, zyffu, znffu, & ! x,y components and norm of the vector: between F points below and above a U point 181 zxffv, zyffv, znffv, & ! x,y components and norm of the vector: between F points left and right a V point 182 zxuuf, zyuuf, znuuf ! x,y components and norm of the vector: between U points below and above a F point 183 !!---------------------------------------------------------------------- 184 185 ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj), & 186 & gsinu(jpi,jpj), gcosu(jpi,jpj), & 187 & gsinv(jpi,jpj), gcosv(jpi,jpj), & 188 & gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr ) 189 IF(lk_mpp) CALL mpp_sum( ierr ) 190 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'angle_msh_geo: unable to allocate arrays' ) 191 192 ! ============================= ! 193 ! Compute the cosinus and sinus ! 194 ! ============================= ! 195 ! (computation done on the north stereographic polar plane) 196 197 DO jj = 2, jpjm1 198 !CDIR NOVERRCHK 199 DO ji = fs_2, jpi ! vector opt. 200 201 ! north pole direction & modulous (at t-point) 202 zlam = glamt(ji,jj) 203 zphi = gphit(ji,jj) 204 zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 205 zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 206 znnpt = zxnpt*zxnpt + zynpt*zynpt 207 208 ! north pole direction & modulous (at u-point) 209 zlam = glamu(ji,jj) 210 zphi = gphiu(ji,jj) 211 zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 212 zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 213 znnpu = zxnpu*zxnpu + zynpu*zynpu 214 215 ! north pole direction & modulous (at v-point) 216 zlam = glamv(ji,jj) 217 zphi = gphiv(ji,jj) 218 zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 219 zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 220 znnpv = zxnpv*zxnpv + zynpv*zynpv 221 222 ! north pole direction & modulous (at f-point) 223 zlam = glamf(ji,jj) 224 zphi = gphif(ji,jj) 225 zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 226 zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 227 znnpf = zxnpf*zxnpf + zynpf*zynpf 228 229 ! j-direction: v-point segment direction (around t-point) 230 zlam = glamv(ji,jj ) 231 zphi = gphiv(ji,jj ) 232 zlan = glamv(ji,jj-1) 233 zphh = gphiv(ji,jj-1) 234 zxvvt = 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) & 235 & - 2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 236 zyvvt = 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) & 237 & - 2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 238 znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt ) ) 239 znvvt = MAX( znvvt, 1.e-14 ) 240 241 ! j-direction: f-point segment direction (around u-point) 242 zlam = glamf(ji,jj ) 243 zphi = gphif(ji,jj ) 244 zlan = glamf(ji,jj-1) 245 zphh = gphif(ji,jj-1) 246 zxffu = 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) & 247 & - 2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 248 zyffu = 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) & 249 & - 2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 250 znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu ) ) 251 znffu = MAX( znffu, 1.e-14 ) 252 253 ! i-direction: f-point segment direction (around v-point) 254 zlam = glamf(ji ,jj) 255 zphi = gphif(ji ,jj) 256 zlan = glamf(ji-1,jj) 257 zphh = gphif(ji-1,jj) 258 zxffv = 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) & 259 & - 2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 260 zyffv = 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) & 261 & - 2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 262 znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv ) ) 263 znffv = MAX( znffv, 1.e-14 ) 264 265 ! j-direction: u-point segment direction (around f-point) 266 zlam = glamu(ji,jj+1) 267 zphi = gphiu(ji,jj+1) 268 zlan = glamu(ji,jj ) 269 zphh = gphiu(ji,jj ) 270 zxuuf = 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) & 271 & - 2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 272 zyuuf = 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) & 273 & - 2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 274 znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf ) ) 275 znuuf = MAX( znuuf, 1.e-14 ) 276 277 ! cosinus and sinus using scalar and vectorial products 278 gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt 279 gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt 280 281 gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu 282 gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu 283 284 gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf 285 gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf 286 287 ! (caution, rotation of 90 degres) 288 gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv 289 gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv 290 291 END DO 292 END DO 293 294 ! =============== ! 295 ! Geographic mesh ! 296 ! =============== ! 297 298 DO jj = 2, jpjm1 299 DO ji = fs_2, jpi ! vector opt. 300 IF( MOD( ABS( glamv(ji,jj) - glamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN 301 gsint(ji,jj) = 0. 302 gcost(ji,jj) = 1. 303 ENDIF 304 IF( MOD( ABS( glamf(ji,jj) - glamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN 305 gsinu(ji,jj) = 0. 306 gcosu(ji,jj) = 1. 307 ENDIF 308 IF( ABS( gphif(ji,jj) - gphif(ji-1,jj) ) < 1.e-8 ) THEN 309 gsinv(ji,jj) = 0. 310 gcosv(ji,jj) = 1. 311 ENDIF 312 IF( MOD( ABS( glamu(ji,jj) - glamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN 313 gsinf(ji,jj) = 0. 314 gcosf(ji,jj) = 1. 315 ENDIF 316 END DO 317 END DO 318 319 ! =========================== ! 320 ! Lateral boundary conditions ! 321 ! =========================== ! 322 323 ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 324 CALL lbc_lnk( gcost, 'T', -1. ) ; CALL lbc_lnk( gsint, 'T', -1. ) 325 CALL lbc_lnk( gcosu, 'U', -1. ) ; CALL lbc_lnk( gsinu, 'U', -1. ) 326 CALL lbc_lnk( gcosv, 'V', -1. ) ; CALL lbc_lnk( gsinv, 'V', -1. ) 327 CALL lbc_lnk( gcosf, 'F', -1. ) ; CALL lbc_lnk( gsinf, 'F', -1. ) 328 329 END SUBROUTINE angle 330 331 332 SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, & 333 pte, ptn ) 334 !!---------------------------------------------------------------------- 335 !! *** ROUTINE geo2oce *** 336 !! 337 !! ** Purpose : 338 !! 339 !! ** Method : Change wind stress from geocentric to east/north 340 !! 341 !! History : 342 !! ! (O. Marti) Original code 343 !! ! 91-03 (G. Madec) 344 !! ! 92-07 (M. Imbard) 345 !! ! 99-11 (M. Imbard) NetCDF format with IOIPSL 346 !! ! 00-08 (D. Ludicone) Reduced section at Bab el Mandeb 347 !! 8.5 ! 02-06 (G. Madec) F90: Free form 348 !! 3.0 ! 07-08 (G. Madec) geo2oce suppress lon/lat agruments 349 !!---------------------------------------------------------------------- 350 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxx, pyy, pzz 351 CHARACTER(len=1) , INTENT(in ) :: cgrid 352 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pte, ptn 353 !! 354 REAL(wp), PARAMETER :: rpi = 3.141592653e0 355 REAL(wp), PARAMETER :: rad = rpi / 180.e0 356 INTEGER :: ig ! 357 !!---------------------------------------------------------------------- 358 359 IF( ALLOCATED( gsinlon ) ) THEN 360 ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , & 361 & gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 362 IF( lk_mpp ) CALL mpp_sum( ierr ) 363 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'angle_msh_geo: unable to allocate arrays' ) 364 ENDIF 365 366 SELECT CASE( cgrid) 367 CASE ( 'T' ) 368 ig = 1 369 IF( .NOT. linit(ig) ) THEN 370 gsinlon(:,:,ig) = SIN( rad * glamt(:,:) ) 371 gcoslon(:,:,ig) = COS( rad * glamt(:,:) ) 372 gsinlat(:,:,ig) = SIN( rad * gphit(:,:) ) 373 gcoslat(:,:,ig) = COS( rad * gphit(:,:) ) 374 linit(ig) = .TRUE. 375 ENDIF 376 CASE ( 'U' ) 377 ig = 2 378 IF( .NOT. linit(ig) ) THEN 379 gsinlon(:,:,ig) = SIN( rad * glamu(:,:) ) 380 gcoslon(:,:,ig) = COS( rad * glamu(:,:) ) 381 gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) ) 382 gcoslat(:,:,ig) = COS( rad * gphiu(:,:) ) 383 linit(ig) = .TRUE. 384 ENDIF 385 CASE ( 'V' ) 386 ig = 3 387 IF( .NOT. linit(ig) ) THEN 388 gsinlon(:,:,ig) = SIN( rad * glamv(:,:) ) 389 gcoslon(:,:,ig) = COS( rad * glamv(:,:) ) 390 gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) ) 391 gcoslat(:,:,ig) = COS( rad * gphiv(:,:) ) 392 linit(ig) = .TRUE. 393 ENDIF 394 CASE ( 'F' ) 395 ig = 4 396 IF( .NOT. linit(ig) ) THEN 397 gsinlon(:,:,ig) = SIN( rad * glamf(:,:) ) 398 gcoslon(:,:,ig) = COS( rad * glamf(:,:) ) 399 gsinlat(:,:,ig) = SIN( rad * gphif(:,:) ) 400 gcoslat(:,:,ig) = COS( rad * gphif(:,:) ) 401 linit(ig) = .TRUE. 402 ENDIF 403 CASE default 404 WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid 405 CALL ctl_stop( ctmp1 ) 406 END SELECT 407 408 pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy 409 ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx & 410 - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy & 411 + gcoslat(:,:,ig) * pzz 412 !!$ ptv = gcoslon(:,:,ig) * gcoslat(:,:,ig) * pxx & 413 !!$ + gsinlon(:,:,ig) * gcoslat(:,:,ig) * pyy & 414 !!$ + gsinlat(:,:,ig) * pzz 358 415 ! 359 END SUBROUTINE repcmo 360 361 362 SUBROUTINE repere( px1, py1, px2, py2, kchoix, cd_type ) 416 END SUBROUTINE geo2oce 417 418 SUBROUTINE oce2geo ( pte, ptn, cgrid, & 419 pxx , pyy , pzz ) 420 !!---------------------------------------------------------------------- 421 !! *** ROUTINE oce2geo *** 422 !! 423 !! ** Purpose : 424 !! 425 !! ** Method : Change vector from east/north to geocentric 426 !! 427 !! History : 428 !! ! (A. Caubel) oce2geo - Original code 429 !!---------------------------------------------------------------------- 430 REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pte, ptn 431 CHARACTER(len=1) , INTENT( IN ) :: cgrid 432 REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ) :: pxx , pyy , pzz 433 !! 434 REAL(wp), PARAMETER :: rpi = 3.141592653E0 435 REAL(wp), PARAMETER :: rad = rpi / 180.e0 436 INTEGER :: ig ! 437 !!---------------------------------------------------------------------- 438 439 IF( 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('STOP', 'angle_msh_geo: 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 493 END SUBROUTINE oce2geo 494 495 496 SUBROUTINE repere ( px1, py1, px2, py2, kchoix, cd_type ) 363 497 !!---------------------------------------------------------------------- 364 498 !! *** ROUTINE repere *** 365 499 !! 366 500 !! ** Purpose : Change vector componantes between a geopgraphic grid 367 !! and a stretched coordinates grid. 368 !!---------------------------------------------------------------------- 369 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: px1, py1 ! the 2 horizontal components to be rotated 370 REAL(wp), DIMENSION(:,:) , INTENT( out) :: px2, py2 ! the 2 horizontal components in the model repere 371 INTEGER , INTENT(in ) :: kchoix ! = 1 change from geographic to model grid. 372 ! ! =-1 change from model to geographic grid 373 CHARACTER(len=1), OPTIONAL, INTENT(in ) :: cd_type ! define the nature of pt2d array grid-points 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 !! * Arguments 515 REAL(wp), INTENT( IN ), DIMENSION(jpi,jpj) :: & 516 px1, py1 ! two horizontal components to be rotated 517 REAL(wp), INTENT( OUT ), DIMENSION(jpi,jpj) :: & 518 px2, py2 ! the two horizontal components in the model repere 519 INTEGER, INTENT( IN ) :: & 520 kchoix ! type of transformation 521 ! = 1 change from geographic to model grid. 522 ! =-1 change from model to geographic grid 523 CHARACTER(len=1), INTENT( IN ), OPTIONAL :: cd_type ! define the nature of pt2d array grid-points 374 524 ! 375 525 CHARACTER(len=1) :: cl_type ! define the nature of pt2d array grid-points (T point by default) 376 526 !!---------------------------------------------------------------------- 377 ! 527 378 528 cl_type = 'T' 379 529 IF( PRESENT(cd_type) ) cl_type = cd_type … … 388 538 CASE DEFAULT ; CALL ctl_stop( 'repere: Syntax Error in the definition of kchoix (1 OR -1' ) 389 539 END SELECT 390 !540 391 541 END SUBROUTINE repere 392 542 393 543 394 SUBROUTINE rot_rep( pxin, pyin, cd_type, cdtodo, prot ) 395 !!---------------------------------------------------------------------- 396 !! *** ROUTINE rot_rep *** 397 !! 398 !! ** Purpose : Rotate the Repere: Change vector componantes between 399 !! geographic grid <--> stretched coordinates grid. 400 !!---------------------------------------------------------------------- 401 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxin, pyin ! vector componantes 402 CHARACTER(len=1), INTENT(in ) :: cd_type ! define the nature of pt2d array grid-points 403 CHARACTER(len=5), INTENT(in ) :: cdtodo ! specify the work to do: 404 !! ! 'en->i' east-north componantes to model i-componante 405 !! ! 'en->j' east-north componantes to model j-componante 406 !! ! 'ij->e' model i-j componantes to east componante 407 !! ! 'ij->n' model i-j componantes to east componante 408 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: prot ! 409 ! 410 INTEGER :: ipt ! 411 !!---------------------------------------------------------------------- 412 413 CALL angle_msh_geo ! initialization of the transformation (just a return if not the 1st call) 414 415 SELECT CASE (cd_type) 416 CASE ('T') ; ipt = 1 417 CASE ('U') ; ipt = 2 418 CASE ('V') ; ipt = 3 419 CASE ('F') ; ipt = 4 420 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 421 END SELECT 422 423 SELECT CASE (cdtodo) 424 CASE ('en->i') ; prot(:,:) = pxin(:,:) * gcos(:,:,ipt) + pyin(:,:) * gsin(:,:,ipt) ! east-north to model i-component 425 CASE ('en->j') ; prot(:,:) = pyin(:,:) * gcos(:,:,ipt) - pxin(:,:) * gsin(:,:,ipt) ! east-north to model j-component 426 CASE ('ij->e') ; prot(:,:) = pxin(:,:) * gcos(:,:,ipt) - pyin(:,:) * gsin(:,:,ipt) ! model i-j to east component 427 CASE ('ij->n') ; prot(:,:) = pyin(:,:) * gcos(:,:,ipt) + pxin(:,:) * gsin(:,:,ipt) ! model i-j to north component 428 CASE DEFAULT ; CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' ) 429 END SELECT 430 431 END SUBROUTINE rot_rep 544 SUBROUTINE obs_rot ( psinu, pcosu, psinv, pcosv ) 545 !!---------------------------------------------------------------------- 546 !! *** ROUTINE obs_rot *** 547 !! 548 !! ** Purpose : Copy gsinu, gcosu, gsinv and gsinv 549 !! to input data for rotations of 550 !! current at observation points 551 !! 552 !! History : 553 !! 9.2 ! 09-02 (K. Mogensen) 554 !!---------------------------------------------------------------------- 555 REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ):: & 556 & psinu, pcosu, psinv, pcosv! copy of data 557 !!---------------------------------------------------------------------- 558 559 ! Initialization of gsin* and gcos* at first call 560 ! ----------------------------------------------- 561 562 IF( lmust_init ) THEN 563 IF(lwp) WRITE(numout,*) 564 IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched' 565 IF(lwp) WRITE(numout,*) ' ~~~~~~~ coordinate transformation' 566 567 CALL angle ! initialization of the transformation 568 lmust_init = .FALSE. 569 570 ENDIF 571 572 psinu(:,:) = gsinu(:,:) 573 pcosu(:,:) = gcosu(:,:) 574 psinv(:,:) = gsinv(:,:) 575 pcosv(:,:) = gcosv(:,:) 576 577 END SUBROUTINE obs_rot 432 578 433 579 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
