- Timestamp:
- 2011-07-11T12:53:56+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcrad.F90
r2715 r2797 5 5 !!================================================================================= 6 6 #if defined key_obc 7 !!---------------------------------------------------------------------------------8 !! obc_rad : call the subroutine for each open boundary9 !! obc_rad_east : compute the east phase velocities10 !! obc_rad_west : compute the west phase velocities11 !! obc_rad_north : compute the north phase velocities12 !! obc_rad_south : compute the south phase velocities13 !!---------------------------------------------------------------------------------14 USE oce ! ocean dynamics and tracers variables15 USE dom_oce ! ocean space and time domain variables16 USE lbclnk ! ocean lateral boundary conditions (or mpp link)17 USE phycst ! physical constants18 USE obc_oce ! ocean open boundary conditions19 USE lib_mpp ! for mppobc20 USE in_out_manager ! I/O units21 22 IMPLICIT NONE23 PRIVATE24 25 PUBLIC obc_rad ! routine called by step.F9026 27 INTEGER :: ji, jj, jk ! dummy loop indices28 29 INTEGER :: & ! ... boundary space indices30 nib = 1, & ! nib = boundary point31 nibm = 2, & ! nibm = 1st interior point32 nibm2 = 3, & ! nibm2 = 2nd interior point33 ! ... boundary time indices34 nit = 1, & ! nit = now35 nitm = 2, & ! nitm = before36 nitm2 = 3 ! nitm2 = before-before37 38 !! * Substitutions39 # include "obc_vectopt_loop_substitute.h90"40 !!---------------------------------------------------------------------------------41 !! NEMO/OPA 3.3 , NEMO Consortium (2010)42 !! $Id$43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)44 !!---------------------------------------------------------------------------------45 46 CONTAINS47 48 SUBROUTINE obc_rad ( kt )49 !!------------------------------------------------------------------------------50 !! SUBROUTINE obc_rad51 !! ********************52 !! ** Purpose :53 !! Perform swap of arrays to calculate radiative phase speeds at the open54 !! boundaries and calculate those phase speeds if the open boundaries are55 !! not fixed. In case of fixed open boundaries does nothing.56 !!57 !! The logical variable lp_obc_east, and/or lp_obc_west, and/or lp_obc_north,58 !! and/or lp_obc_south allow the user to determine which boundary is an59 !! open one (must be done in the param_obc.h90 file).60 !!61 !! ** Reference :62 !! Marchesiello P., 1995, these de l'universite J. Fourier, Grenoble, France.63 !!64 !! History :65 !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90 from the66 !! J. Molines and G. Madec version67 !!------------------------------------------------------------------------------68 INTEGER, INTENT( in ) :: kt69 !!----------------------------------------------------------------------70 71 IF( lp_obc_east .AND. .NOT.lfbceast ) CALL obc_rad_east ( kt ) ! East open boundary72 73 IF( lp_obc_west .AND. .NOT.lfbcwest ) CALL obc_rad_west ( kt ) ! West open boundary74 75 IF( lp_obc_north .AND. .NOT.lfbcnorth ) CALL obc_rad_north( kt ) ! North open boundary76 77 IF( lp_obc_south .AND. .NOT.lfbcsouth ) CALL obc_rad_south( kt ) ! South open boundary78 79 END SUBROUTINE obc_rad80 81 82 SUBROUTINE obc_rad_east ( kt )83 !!------------------------------------------------------------------------------84 !! *** SUBROUTINE obc_rad_east ***85 !!86 !! ** Purpose :87 !! Perform swap of arrays to calculate radiative phase speeds at the open88 !! east boundary and calculate those phase speeds if this OBC is not fixed.89 !! In case of fixed OBC, this subrountine is not called.90 !!91 !! History :92 !! ! 95-03 (J.-M. Molines) Original from SPEM93 !! ! 97-07 (G. Madec, J.-M. Molines) additions94 !! ! 97-12 (M. Imbard) Mpp adaptation95 !! ! 00-06 (J.-M. Molines)96 !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F9097 !!------------------------------------------------------------------------------98 !! * Arguments99 INTEGER, INTENT( in ) :: kt100 101 !! * Local declarations102 INTEGER :: ij103 REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy104 REAL(wp) :: zucb, zucbm, zucbm2105 !!------------------------------------------------------------------------------106 107 ! 1. Swap arrays before calculating radiative velocities108 ! ------------------------------------------------------109 110 ! 1.1 zonal velocity111 ! -------------------112 113 IF( kt > nit000 .OR. ln_rstart ) THEN114 115 ! ... advance in time (time filter, array swap)116 DO jk = 1, jpkm1117 DO jj = 1, jpj118 uebnd(jj,jk,nib ,nitm2) = uebnd(jj,jk,nib ,nitm)*uemsk(jj,jk)119 uebnd(jj,jk,nibm ,nitm2) = uebnd(jj,jk,nibm ,nitm)*uemsk(jj,jk)120 uebnd(jj,jk,nibm2,nitm2) = uebnd(jj,jk,nibm2,nitm)*uemsk(jj,jk)121 END DO122 END DO123 ! ... fields nitm <== nit plus time filter at the boundary124 DO ji = fs_nie0, fs_nie1 ! Vector opt.125 DO jk = 1, jpkm1126 DO jj = 1, jpj127 uebnd(jj,jk,nib ,nitm) = uebnd(jj,jk,nib, nit)*uemsk(jj,jk)128 uebnd(jj,jk,nibm ,nitm) = uebnd(jj,jk,nibm ,nit)*uemsk(jj,jk)129 uebnd(jj,jk,nibm2,nitm) = uebnd(jj,jk,nibm2,nit)*uemsk(jj,jk)130 ! ... fields nit <== now (kt+1)131 ! ... Total or baroclinic velocity at b, bm and bm2132 zucb = un(ji,jj,jk)133 zucbm = un(ji-1,jj,jk)134 zucbm2 = un(ji-2,jj,jk)135 uebnd(jj,jk,nib ,nit) = zucb *uemsk(jj,jk)136 uebnd(jj,jk,nibm ,nit) = zucbm *uemsk(jj,jk)137 uebnd(jj,jk,nibm2,nit) = zucbm2 *uemsk(jj,jk)138 END DO139 END DO140 END DO141 IF( lk_mpp ) CALL mppobc(uebnd,jpjed,jpjef,jpieob,jpk*3*3,2,jpj, numout )142 143 ! ... extremeties nie0, nie1144 ij = jpjed +1 - njmpp145 IF( ij >= 2 .AND. ij < jpjm1 ) THEN146 DO jk = 1,jpkm1147 uebnd(ij,jk,nibm,nitm) = uebnd(ij+1 ,jk,nibm,nitm)148 END DO149 END IF150 ij = jpjef +1 - njmpp151 IF( ij >= 2 .AND. ij < jpjm1 ) THEN152 DO jk = 1,jpkm1153 uebnd(ij,jk,nibm,nitm) = uebnd(ij-1 ,jk,nibm,nitm)154 END DO155 END IF156 157 ! 1.2 tangential velocity158 ! -----------------------159 160 ! ... advance in time (time filter, array swap)161 DO jk = 1, jpkm1162 DO jj = 1, jpj163 ! ... fields nitm2 <== nitm164 vebnd(jj,jk,nib ,nitm2) = vebnd(jj,jk,nib ,nitm)*vemsk(jj,jk)165 vebnd(jj,jk,nibm ,nitm2) = vebnd(jj,jk,nibm ,nitm)*vemsk(jj,jk)166 vebnd(jj,jk,nibm2,nitm2) = vebnd(jj,jk,nibm2,nitm)*vemsk(jj,jk)167 END DO168 END DO169 170 DO ji = fs_nie0+1, fs_nie1+1 ! Vector opt.171 DO jk = 1, jpkm1172 DO jj = 1, jpj173 vebnd(jj,jk,nib ,nitm) = vebnd(jj,jk,nib, nit)*vemsk(jj,jk)174 vebnd(jj,jk,nibm ,nitm) = vebnd(jj,jk,nibm ,nit)*vemsk(jj,jk)175 vebnd(jj,jk,nibm2,nitm) = vebnd(jj,jk,nibm2,nit)*vemsk(jj,jk)176 ! ... fields nit <== now (kt+1)177 vebnd(jj,jk,nib ,nit) = vn(ji ,jj,jk)*vemsk(jj,jk)178 vebnd(jj,jk,nibm ,nit) = vn(ji-1,jj,jk)*vemsk(jj,jk)179 vebnd(jj,jk,nibm2,nit) = vn(ji-2,jj,jk)*vemsk(jj,jk)180 END DO181 END DO182 END DO183 IF( lk_mpp ) CALL mppobc(vebnd,jpjed,jpjef,jpieob+1,jpk*3*3,2,jpj, numout )184 185 !... extremeties nie0, nie1186 ij = jpjed +1 - njmpp187 IF( ij >= 2 .AND. ij < jpjm1 ) THEN188 DO jk = 1,jpkm1189 vebnd(ij,jk,nibm,nitm) = vebnd(ij+1 ,jk,nibm,nitm)190 END DO191 END IF192 ij = jpjef +1 - njmpp193 IF( ij >= 2 .AND. ij < jpjm1 ) THEN194 DO jk = 1,jpkm1195 vebnd(ij,jk,nibm,nitm) = vebnd(ij-1 ,jk,nibm,nitm)196 END DO197 END IF198 199 ! 1.3 Temperature and salinity200 ! ----------------------------201 202 ! ... advance in time (time filter, array swap)203 DO jk = 1, jpkm1204 DO jj = 1, jpj205 ! ... fields nitm <== nit plus time filter at the boundary206 tebnd(jj,jk,nib,nitm) = tebnd(jj,jk,nib,nit)*temsk(jj,jk)207 sebnd(jj,jk,nib,nitm) = sebnd(jj,jk,nib,nit)*temsk(jj,jk)208 END DO209 END DO210 211 DO ji = fs_nie0+1, fs_nie1+1 ! Vector opt.212 DO jk = 1, jpkm1213 DO jj = 1, jpj214 tebnd(jj,jk,nibm,nitm) = tebnd(jj,jk,nibm,nit)*temsk(jj,jk)215 sebnd(jj,jk,nibm,nitm) = sebnd(jj,jk,nibm,nit)*temsk(jj,jk)216 ! ... fields nit <== now (kt+1)217 tebnd(jj,jk,nib ,nit) = tn(ji ,jj,jk)*temsk(jj,jk)218 tebnd(jj,jk,nibm ,nit) = tn(ji-1,jj,jk)*temsk(jj,jk)219 sebnd(jj,jk,nib ,nit) = sn(ji ,jj,jk)*temsk(jj,jk)220 sebnd(jj,jk,nibm ,nit) = sn(ji-1,jj,jk)*temsk(jj,jk)221 END DO222 END DO223 END DO224 IF( lk_mpp ) CALL mppobc(tebnd,jpjed,jpjef,jpieob+1,jpk*2*2,2,jpj, numout )225 IF( lk_mpp ) CALL mppobc(sebnd,jpjed,jpjef,jpieob+1,jpk*2*2,2,jpj, numout )226 227 ! ... extremeties nie0, nie1228 ij = jpjed +1 - njmpp229 IF( ij >= 2 .AND. ij < jpjm1 ) THEN230 DO jk = 1,jpkm1231 tebnd(ij,jk,nibm,nitm) = tebnd(ij+1 ,jk,nibm,nitm)232 sebnd(ij,jk,nibm,nitm) = sebnd(ij+1 ,jk,nibm,nitm)233 END DO234 END IF235 ij = jpjef +1 - njmpp236 IF( ij >= 2 .AND. ij < jpjm1 ) THEN237 DO jk = 1,jpkm1238 tebnd(ij,jk,nibm,nitm) = tebnd(ij-1 ,jk,nibm,nitm)239 sebnd(ij,jk,nibm,nitm) = sebnd(ij-1 ,jk,nibm,nitm)240 END DO241 END IF242 243 END IF ! End of array swap244 245 ! 2 - Calculation of radiation velocities246 ! ---------------------------------------247 248 IF( kt >= nit000 +3 .OR. ln_rstart ) THEN249 250 ! 2.1 Calculate the normal velocity U based on phase velocity u_cxebnd251 ! ---------------------------------------------------------------------252 !253 ! nibm2 nibm nib254 ! | nibm | nib |///255 ! | | | | |///256 ! jj-line --f----v----f----v----f---257 ! | | | | |///258 ! | | |///259 ! jj-line u T u T u///260 ! | | |///261 ! | | | | |///262 ! jpieob-2 jpieob-1 jpieob263 ! | |264 ! jpieob-1 jpieob265 !266 ! ... (jpjedp1, jpjefm1),jpieob267 DO ji = fs_nie0, fs_nie1 ! Vector opt.268 DO jk = 1, jpkm1269 DO jj = 2, jpjm1270 ! ... 2* gradi(u) (T-point i=nibm, time mean)271 z2dx = ( uebnd(jj,jk,nibm ,nit) + uebnd(jj,jk,nibm ,nitm2) &272 - 2.*uebnd(jj,jk,nibm2,nitm) ) / e1t(ji-1,jj)273 ! ... 2* gradj(u) (u-point i=nibm, time nitm)274 z2dy = ( uebnd(jj+1,jk,nibm,nitm) - uebnd(jj-1,jk,nibm,nitm) ) / e2u(ji-1,jj)275 ! ... square of the norm of grad(u)276 z4nor2 = z2dx * z2dx + z2dy * z2dy277 ! ... minus time derivative (leap-frog) at nibm, without / 2 dt278 zdt = uebnd(jj,jk,nibm,nitm2) - uebnd(jj,jk,nibm,nit)279 ! ... i-phase speed ratio (bounded by 1)280 IF( z4nor2 == 0. ) THEN281 z4nor2=.00001282 END IF283 z05cx = zdt * z2dx / z4nor2284 u_cxebnd(jj,jk) = z05cx*uemsk(jj,jk)285 END DO286 END DO287 END DO288 289 ! 2.2 Calculate the tangential velocity based on phase velocity v_cxebnd290 ! -----------------------------------------------------------------------291 !292 ! nibm2 nibm nib293 ! | nibm | nib///|///294 ! | | | |////|///295 ! jj-line --v----f----v----f----v---296 ! | | | |////|///297 ! | | | |////|///298 ! | jpieob-1| jpieob /|///299 ! | | |300 ! jpieob-1 jpieob jpieob+1301 !302 ! ... (jpjedp1, jpjefm1), jpieob+1303 DO ji = fs_nie0+1, fs_nie1+1 ! Vector opt.304 DO jk = 1, jpkm1305 DO jj = 2, jpjm1306 ! ... 2* i-gradient of v (f-point i=nibm, time mean)307 z2dx = ( vebnd(jj,jk,nibm ,nit) + vebnd(jj,jk,nibm ,nitm2) &308 - 2.*vebnd(jj,jk,nibm2,nitm) ) / e1f(ji-2,jj)309 ! ... 2* j-gradient of v (v-point i=nibm, time nitm)310 z2dy = ( vebnd(jj+1,jk,nibm,nitm) - vebnd(jj-1,jk,nibm,nitm) ) / e2v(ji-1,jj)311 ! ... square of the norm of grad(v)312 z4nor2 = z2dx * z2dx + z2dy * z2dy313 ! ... minus time derivative (leap-frog) at nibm, without / 2 dt314 zdt = vebnd(jj,jk,nibm,nitm2) - vebnd(jj,jk,nibm,nit)315 ! ... i-phase speed ratio (bounded by 1) and save the unbounded phase316 ! velocity ratio no divided by e1f for the tracer radiation317 IF( z4nor2 == 0. ) THEN318 z4nor2=.000001319 END IF320 z05cx = zdt * z2dx / z4nor2321 v_cxebnd(jj,jk) = z05cx*vemsk(jj,jk)322 END DO323 END DO324 END DO325 IF( lk_mpp ) CALL mppobc(v_cxebnd,jpjed,jpjef,jpieob+1,jpk,2,jpj, numout )326 327 ! ... extremeties nie0, nie1328 ij = jpjed +1 - njmpp329 IF( ij >= 2 .AND. ij < jpjm1 ) THEN330 DO jk = 1,jpkm1331 v_cxebnd(ij,jk) = v_cxebnd(ij+1 ,jk)332 END DO333 END IF334 ij = jpjef +1 - njmpp335 IF( ij >= 2 .AND. ij < jpjm1 ) THEN336 DO jk = 1,jpkm1337 v_cxebnd(ij,jk) = v_cxebnd(ij-1 ,jk)338 END DO339 END IF340 341 END IF342 343 END SUBROUTINE obc_rad_east344 345 346 SUBROUTINE obc_rad_west ( kt )347 !!------------------------------------------------------------------------------348 !! *** SUBROUTINE obc_rad_west ***349 !!350 !! ** Purpose :351 !! Perform swap of arrays to calculate radiative phase speeds at the open352 !! west boundary and calculate those phase speeds if this OBC is not fixed.353 !! In case of fixed OBC, this subrountine is not called.354 !!355 !! History :356 !! ! 95-03 (J.-M. Molines) Original from SPEM357 !! ! 97-07 (G. Madec, J.-M. Molines) additions358 !! ! 97-12 (M. Imbard) Mpp adaptation359 !! ! 00-06 (J.-M. Molines)360 !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90361 !!------------------------------------------------------------------------------362 !! * Arguments363 INTEGER, INTENT( in ) :: kt364 365 !! * Local declarations366 INTEGER :: ij367 REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy368 REAL(wp) :: zucb, zucbm, zucbm2369 !!------------------------------------------------------------------------------370 371 ! 1. Swap arrays before calculating radiative velocities372 ! ------------------------------------------------------373 374 ! 1.1 zonal velocity375 ! -------------------376 377 IF( kt > nit000 .OR. ln_rstart ) THEN378 379 ! ... advance in time (time filter, array swap)380 DO jk = 1, jpkm1381 DO jj = 1, jpj382 uwbnd(jj,jk,nib ,nitm2) = uwbnd(jj,jk,nib ,nitm)*uwmsk(jj,jk)383 uwbnd(jj,jk,nibm ,nitm2) = uwbnd(jj,jk,nibm ,nitm)*uwmsk(jj,jk)384 uwbnd(jj,jk,nibm2,nitm2) = uwbnd(jj,jk,nibm2,nitm)*uwmsk(jj,jk)385 END DO386 END DO387 388 ! ... fields nitm <== nit plus time filter at the boundary389 DO ji = fs_niw0, fs_niw1 ! Vector opt.390 DO jk = 1, jpkm1391 DO jj = 1, jpj392 uwbnd(jj,jk,nib ,nitm) = uwbnd(jj,jk,nib ,nit)*uwmsk(jj,jk)393 uwbnd(jj,jk,nibm ,nitm) = uwbnd(jj,jk,nibm ,nit)*uwmsk(jj,jk)394 uwbnd(jj,jk,nibm2,nitm) = uwbnd(jj,jk,nibm2,nit)*uwmsk(jj,jk)395 ! ... total or baroclinic velocity at b, bm and bm2396 zucb = un (ji,jj,jk)397 zucbm = un (ji+1,jj,jk)398 zucbm2 = un (ji+2,jj,jk)399 400 ! ... fields nit <== now (kt+1)401 uwbnd(jj,jk,nib ,nit) = zucb *uwmsk(jj,jk)402 uwbnd(jj,jk,nibm ,nit) = zucbm *uwmsk(jj,jk)403 uwbnd(jj,jk,nibm2,nit) = zucbm2*uwmsk(jj,jk)404 END DO405 END DO406 END DO407 IF( lk_mpp ) CALL mppobc(uwbnd,jpjwd,jpjwf,jpiwob,jpk*3*3,2,jpj, numout )408 409 ! ... extremeties niw0, niw1410 ij = jpjwd +1 - njmpp411 IF( ij >= 2 .AND. ij < jpjm1 ) THEN412 DO jk = 1,jpkm1413 uwbnd(ij,jk,nibm,nitm) = uwbnd(ij+1 ,jk,nibm,nitm)414 END DO415 END IF416 ij = jpjwf +1 - njmpp417 IF( ij >= 2 .AND. ij < jpjm1 ) THEN418 DO jk = 1,jpkm1419 uwbnd(ij,jk,nibm,nitm) = uwbnd(ij-1 ,jk,nibm,nitm)420 END DO421 END IF422 423 ! 1.2 tangential velocity424 ! -----------------------425 426 ! ... advance in time (time filter, array swap)427 DO jk = 1, jpkm1428 DO jj = 1, jpj429 ! ... fields nitm2 <== nitm430 vwbnd(jj,jk,nib ,nitm2) = vwbnd(jj,jk,nib ,nitm)*vwmsk(jj,jk)431 vwbnd(jj,jk,nibm ,nitm2) = vwbnd(jj,jk,nibm ,nitm)*vwmsk(jj,jk)432 vwbnd(jj,jk,nibm2,nitm2) = vwbnd(jj,jk,nibm2,nitm)*vwmsk(jj,jk)433 END DO434 END DO435 436 DO ji = fs_niw0, fs_niw1 ! Vector opt.437 DO jk = 1, jpkm1438 DO jj = 1, jpj439 vwbnd(jj,jk,nib ,nitm) = vwbnd(jj,jk,nib, nit)*vwmsk(jj,jk)440 vwbnd(jj,jk,nibm ,nitm) = vwbnd(jj,jk,nibm ,nit)*vwmsk(jj,jk)441 vwbnd(jj,jk,nibm2,nitm) = vwbnd(jj,jk,nibm2,nit)*vwmsk(jj,jk)442 ! ... fields nit <== now (kt+1)443 vwbnd(jj,jk,nib ,nit) = vn(ji ,jj,jk)*vwmsk(jj,jk)444 vwbnd(jj,jk,nibm ,nit) = vn(ji+1,jj,jk)*vwmsk(jj,jk)445 vwbnd(jj,jk,nibm2,nit) = vn(ji+2,jj,jk)*vwmsk(jj,jk)446 END DO447 END DO448 END DO449 IF( lk_mpp ) CALL mppobc(vwbnd,jpjwd,jpjwf,jpiwob,jpk*3*3,2,jpj, numout )450 451 ! ... extremeties niw0, niw1452 ij = jpjwd +1 - njmpp453 IF( ij >= 2 .AND. ij < jpjm1 ) THEN454 DO jk = 1,jpkm1455 vwbnd(ij,jk,nibm,nitm) = vwbnd(ij+1 ,jk,nibm,nitm)456 END DO457 END IF458 ij = jpjwf +1 - njmpp459 IF( ij >= 2 .AND. ij < jpjm1 ) THEN460 DO jk = 1,jpkm1461 vwbnd(ij,jk,nibm,nitm) = vwbnd(ij-1 ,jk,nibm,nitm)462 END DO463 END IF464 465 ! 1.3 Temperature and salinity466 ! ----------------------------467 468 ! ... advance in time (time filter, array swap)469 DO jk = 1, jpkm1470 DO jj = 1, jpj471 ! ... fields nitm <== nit plus time filter at the boundary472 twbnd(jj,jk,nib,nitm) = twbnd(jj,jk,nib,nit)*twmsk(jj,jk)473 swbnd(jj,jk,nib,nitm) = swbnd(jj,jk,nib,nit)*twmsk(jj,jk)474 END DO475 END DO476 477 DO ji = fs_niw0, fs_niw1 ! Vector opt.478 DO jk = 1, jpkm1479 DO jj = 1, jpj480 twbnd(jj,jk,nibm ,nitm) = twbnd(jj,jk,nibm ,nit)*twmsk(jj,jk)481 swbnd(jj,jk,nibm ,nitm) = swbnd(jj,jk,nibm ,nit)*twmsk(jj,jk)482 ! ... fields nit <== now (kt+1)483 twbnd(jj,jk,nib ,nit) = tn(ji ,jj,jk)*twmsk(jj,jk)484 twbnd(jj,jk,nibm ,nit) = tn(ji+1 ,jj,jk)*twmsk(jj,jk)485 swbnd(jj,jk,nib ,nit) = sn(ji ,jj,jk)*twmsk(jj,jk)486 swbnd(jj,jk,nibm ,nit) = sn(ji+1 ,jj,jk)*twmsk(jj,jk)487 END DO488 END DO489 END DO490 IF( lk_mpp ) CALL mppobc(twbnd,jpjwd,jpjwf,jpiwob,jpk*2*2,2,jpj, numout )491 IF( lk_mpp ) CALL mppobc(swbnd,jpjwd,jpjwf,jpiwob,jpk*2*2,2,jpj, numout )492 493 ! ... extremeties niw0, niw1494 ij = jpjwd +1 - njmpp495 IF( ij >= 2 .AND. ij < jpjm1 ) THEN496 DO jk = 1,jpkm1497 twbnd(ij,jk,nibm,nitm) = twbnd(ij+1 ,jk,nibm,nitm)498 swbnd(ij,jk,nibm,nitm) = swbnd(ij+1 ,jk,nibm,nitm)499 END DO500 END IF501 ij = jpjwf +1 - njmpp502 IF( ij >= 2 .AND. ij < jpjm1 ) THEN503 DO jk = 1,jpkm1504 twbnd(ij,jk,nibm,nitm) = twbnd(ij-1 ,jk,nibm,nitm)505 swbnd(ij,jk,nibm,nitm) = swbnd(ij-1 ,jk,nibm,nitm)506 END DO507 END IF508 509 END IF ! End of array swap510 511 ! 2 - Calculation of radiation velocities512 ! ---------------------------------------513 514 IF( kt >= nit000 +3 .OR. ln_rstart ) THEN515 516 ! 2.1 Calculate the normal velocity U based on phase velocity u_cxwbnd517 ! ----------------------------------------------------------------------518 !519 ! nib nibm nibm2520 ! ///| nib | nibm |521 ! ///| | | | |522 ! ---f----v----f----v----f-- jj-line523 ! ///| | | | |524 ! ///| | |525 ! ///u T u T u jj-line526 ! ///| | |527 ! ///| | | | |528 ! jpiwob jpiwob+1 jpiwob+2529 ! | |530 ! jpiwob+1 jpiwob+2531 !532 ! ... If free surface formulation:533 ! ... radiative conditions on the total part + relaxation toward climatology534 ! ... (jpjwdp1, jpjwfm1), jpiwob535 DO ji = fs_niw0, fs_niw1 ! Vector opt.536 DO jk = 1, jpkm1537 DO jj = 2, jpjm1538 ! ... 2* gradi(u) (T-point i=nibm, time mean)539 z2dx = ( - uwbnd(jj,jk,nibm ,nit) - uwbnd(jj,jk,nibm ,nitm2) &540 + 2.*uwbnd(jj,jk,nibm2,nitm) ) / e1t(ji+2,jj)541 ! ... 2* gradj(u) (u-point i=nibm, time nitm)542 z2dy = ( uwbnd(jj+1,jk,nibm,nitm) - uwbnd(jj-1,jk,nibm,nitm) ) / e2u(ji+1,jj)543 ! ... square of the norm of grad(u)544 z4nor2 = z2dx * z2dx + z2dy * z2dy545 ! ... minus time derivative (leap-frog) at nibm, without / 2 dt546 zdt = uwbnd(jj,jk,nibm,nitm2) - uwbnd(jj,jk,nibm,nit)547 ! ... i-phase speed ratio (bounded by -1)548 IF( z4nor2 == 0. ) THEN549 z4nor2=0.00001550 END IF551 z05cx = zdt * z2dx / z4nor2552 u_cxwbnd(jj,jk)=z05cx*uwmsk(jj,jk)553 END DO554 END DO555 END DO556 557 ! 2.2 Calculate the tangential velocity based on phase velocity v_cxwbnd558 ! -----------------------------------------------------------------------559 !560 ! nib nibm nibm2561 ! ///|///nib | nibm | nibm2562 ! ///|////| | | | | |563 ! ---v----f----v----f----v----f----v-- jj-line564 ! ///|////| | | | | |565 ! ///|////| | | | | |566 ! jpiwob jpiwob+1 jpiwob+2567 ! | | |568 ! jpiwob jpiwob+1 jpiwob+2569 !570 ! ... radiative condition plus Raymond-Kuo571 ! ... (jpjwdp1, jpjwfm1),jpiwob572 DO ji = fs_niw0, fs_niw1 ! Vector opt.573 DO jk = 1, jpkm1574 DO jj = 2, jpjm1575 ! ... 2* i-gradient of v (f-point i=nibm, time mean)576 z2dx = ( - vwbnd(jj,jk,nibm ,nit) - vwbnd(jj,jk,nibm ,nitm2) &577 + 2.*vwbnd(jj,jk,nibm2,nitm) ) / e1f(ji+1,jj)578 ! ... 2* j-gradient of v (v-point i=nibm, time nitm)579 z2dy = ( vwbnd(jj+1,jk,nibm,nitm) - vwbnd(jj-1,jk,nibm,nitm) ) / e2v(ji+1,jj)580 ! ... square of the norm of grad(v)581 z4nor2 = z2dx * z2dx + z2dy * z2dy582 ! ... minus time derivative (leap-frog) at nibm, without / 2 dt583 zdt = vwbnd(jj,jk,nibm,nitm2) - vwbnd(jj,jk,nibm,nit)584 ! ... i-phase speed ratio (bounded by -1) and save the unbounded phase585 ! velocity ratio no divided by e1f for the tracer radiation586 IF( z4nor2 == 0) THEN587 z4nor2=0.000001588 endif589 z05cx = zdt * z2dx / z4nor2590 v_cxwbnd(jj,jk) = z05cx*vwmsk(jj,jk)591 END DO592 END DO593 END DO594 IF( lk_mpp ) CALL mppobc(v_cxwbnd,jpjwd,jpjwf,jpiwob,jpk,2,jpj, numout )595 596 ! ... extremeties niw0, niw1597 ij = jpjwd +1 - njmpp598 IF( ij >= 2 .AND. ij < jpjm1 ) THEN599 DO jk = 1,jpkm1600 v_cxwbnd(ij,jk) = v_cxwbnd(ij+1 ,jk)601 END DO602 END IF603 ij = jpjwf +1 - njmpp604 IF( ij >= 2 .AND. ij < jpjm1 ) THEN605 DO jk = 1,jpkm1606 v_cxwbnd(ij,jk) = v_cxwbnd(ij-1 ,jk)607 END DO608 END IF609 610 END IF611 612 END SUBROUTINE obc_rad_west613 614 615 SUBROUTINE obc_rad_north ( kt )616 !!------------------------------------------------------------------------------617 !! *** SUBROUTINE obc_rad_north ***618 !!619 !! ** Purpose :620 !! Perform swap of arrays to calculate radiative phase speeds at the open621 !! north boundary and calculate those phase speeds if this OBC is not fixed.622 !! In case of fixed OBC, this subrountine is not called.623 !!624 !! History :625 !! ! 95-03 (J.-M. Molines) Original from SPEM626 !! ! 97-07 (G. Madec, J.-M. Molines) additions627 !! ! 97-12 (M. Imbard) Mpp adaptation628 !! ! 00-06 (J.-M. Molines)629 !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90630 !!------------------------------------------------------------------------------631 !! * Arguments632 INTEGER, INTENT( in ) :: kt633 634 !! * Local declarations635 INTEGER :: ii636 REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy637 REAL(wp) :: zvcb, zvcbm, zvcbm2638 !!------------------------------------------------------------------------------639 640 ! 1. Swap arrays before calculating radiative velocities641 ! ------------------------------------------------------642 643 ! 1.1 zonal velocity644 ! -------------------645 646 IF( kt > nit000 .OR. ln_rstart ) THEN647 648 ! ... advance in time (time filter, array swap)649 DO jk = 1, jpkm1650 DO ji = 1, jpi651 ! ... fields nitm2 <== nitm652 unbnd(ji,jk,nib ,nitm2) = unbnd(ji,jk,nib ,nitm)*unmsk(ji,jk)653 unbnd(ji,jk,nibm ,nitm2) = unbnd(ji,jk,nibm ,nitm)*unmsk(ji,jk)654 unbnd(ji,jk,nibm2,nitm2) = unbnd(ji,jk,nibm2,nitm)*unmsk(ji,jk)655 END DO656 END DO657 658 DO jj = fs_njn0+1, fs_njn1+1 ! Vector opt.659 DO jk = 1, jpkm1660 DO ji = 1, jpi661 unbnd(ji,jk,nib ,nitm) = unbnd(ji,jk,nib, nit)*unmsk(ji,jk)662 unbnd(ji,jk,nibm ,nitm) = unbnd(ji,jk,nibm ,nit)*unmsk(ji,jk)663 unbnd(ji,jk,nibm2,nitm) = unbnd(ji,jk,nibm2,nit)*unmsk(ji,jk)664 ! ... fields nit <== now (kt+1)665 unbnd(ji,jk,nib ,nit) = un(ji,jj, jk)*unmsk(ji,jk)666 unbnd(ji,jk,nibm ,nit) = un(ji,jj-1,jk)*unmsk(ji,jk)667 unbnd(ji,jk,nibm2,nit) = un(ji,jj-2,jk)*unmsk(ji,jk)668 END DO669 END DO670 END DO671 IF( lk_mpp ) CALL mppobc(unbnd,jpind,jpinf,jpjnob+1,jpk*3*3,1,jpi, numout )672 673 ! ... extremeties njn0,njn1674 ii = jpind + 1 - nimpp675 IF( ii >= 2 .AND. ii < jpim1 ) THEN676 DO jk = 1, jpkm1677 unbnd(ii,jk,nibm,nitm) = unbnd(ii+1,jk,nibm,nitm)678 END DO679 END IF680 ii = jpinf + 1 - nimpp681 IF( ii >= 2 .AND. ii < jpim1 ) THEN682 DO jk = 1, jpkm1683 unbnd(ii,jk,nibm,nitm) = unbnd(ii-1,jk,nibm,nitm)684 END DO685 END IF686 687 ! 1.2. normal velocity688 ! --------------------689 690 ! ... advance in time (time filter, array swap)691 DO jk = 1, jpkm1692 DO ji = 1, jpi693 ! ... fields nitm2 <== nitm694 vnbnd(ji,jk,nib ,nitm2) = vnbnd(ji,jk,nib ,nitm)*vnmsk(ji,jk)695 vnbnd(ji,jk,nibm ,nitm2) = vnbnd(ji,jk,nibm ,nitm)*vnmsk(ji,jk)696 vnbnd(ji,jk,nibm2,nitm2) = vnbnd(ji,jk,nibm2,nitm)*vnmsk(ji,jk)697 END DO698 END DO699 700 DO jj = fs_njn0, fs_njn1 ! Vector opt.701 DO jk = 1, jpkm1702 DO ji = 1, jpi703 vnbnd(ji,jk,nib ,nitm) = vnbnd(ji,jk,nib, nit)*vnmsk(ji,jk)704 vnbnd(ji,jk,nibm ,nitm) = vnbnd(ji,jk,nibm ,nit)*vnmsk(ji,jk)705 vnbnd(ji,jk,nibm2,nitm) = vnbnd(ji,jk,nibm2,nit)*vnmsk(ji,jk)706 ! ... fields nit <== now (kt+1)707 ! ... total or baroclinic velocity at b, bm and bm2708 zvcb = vn (ji,jj,jk)709 zvcbm = vn (ji,jj-1,jk)710 zvcbm2 = vn (ji,jj-2,jk)711 ! ... fields nit <== now (kt+1)712 vnbnd(ji,jk,nib ,nit) = zvcb *vnmsk(ji,jk)713 vnbnd(ji,jk,nibm ,nit) = zvcbm *vnmsk(ji,jk)714 vnbnd(ji,jk,nibm2,nit) = zvcbm2*vnmsk(ji,jk)715 END DO716 END DO717 END DO718 IF( lk_mpp ) CALL mppobc(vnbnd,jpind,jpinf,jpjnob,jpk*3*3,1,jpi, numout )719 720 ! ... extremeties njn0,njn1721 ii = jpind + 1 - nimpp722 IF( ii >= 2 .AND. ii < jpim1 ) THEN723 DO jk = 1, jpkm1724 vnbnd(ii,jk,nibm,nitm) = vnbnd(ii+1,jk,nibm,nitm)725 END DO726 END IF727 ii = jpinf + 1 - nimpp728 IF( ii >= 2 .AND. ii < jpim1 ) THEN729 DO jk = 1, jpkm1730 vnbnd(ii,jk,nibm,nitm) = vnbnd(ii-1,jk,nibm,nitm)731 END DO732 END IF733 734 ! 1.3 Temperature and salinity735 ! ----------------------------736 737 ! ... advance in time (time filter, array swap)738 DO jk = 1, jpkm1739 DO ji = 1, jpi740 ! ... fields nitm <== nit plus time filter at the boundary741 tnbnd(ji,jk,nib ,nitm) = tnbnd(ji,jk,nib,nit)*tnmsk(ji,jk)742 snbnd(ji,jk,nib ,nitm) = snbnd(ji,jk,nib,nit)*tnmsk(ji,jk)743 END DO744 END DO745 746 DO jj = fs_njn0+1, fs_njn1+1 ! Vector opt.747 DO jk = 1, jpkm1748 DO ji = 1, jpi749 tnbnd(ji,jk,nibm ,nitm) = tnbnd(ji,jk,nibm ,nit)*tnmsk(ji,jk)750 snbnd(ji,jk,nibm ,nitm) = snbnd(ji,jk,nibm ,nit)*tnmsk(ji,jk)751 ! ... fields nit <== now (kt+1)752 tnbnd(ji,jk,nib ,nit) = tn(ji,jj, jk)*tnmsk(ji,jk)753 tnbnd(ji,jk,nibm ,nit) = tn(ji,jj-1,jk)*tnmsk(ji,jk)754 snbnd(ji,jk,nib ,nit) = sn(ji,jj, jk)*tnmsk(ji,jk)755 snbnd(ji,jk,nibm ,nit) = sn(ji,jj-1,jk)*tnmsk(ji,jk)756 END DO757 END DO758 END DO759 IF( lk_mpp ) CALL mppobc(tnbnd,jpind,jpinf,jpjnob+1,jpk*2*2,1,jpi, numout )760 IF( lk_mpp ) CALL mppobc(snbnd,jpind,jpinf,jpjnob+1,jpk*2*2,1,jpi, numout )761 762 ! ... extremeties njn0,njn1763 ii = jpind + 1 - nimpp764 IF( ii >= 2 .AND. ii < jpim1 ) THEN765 DO jk = 1, jpkm1766 tnbnd(ii,jk,nibm,nitm) = tnbnd(ii+1,jk,nibm,nitm)767 snbnd(ii,jk,nibm,nitm) = snbnd(ii+1,jk,nibm,nitm)768 END DO769 END IF770 ii = jpinf + 1 - nimpp771 IF( ii >= 2 .AND. ii < jpim1 ) THEN772 DO jk = 1, jpkm1773 tnbnd(ii,jk,nibm,nitm) = tnbnd(ii-1,jk,nibm,nitm)774 snbnd(ii,jk,nibm,nitm) = snbnd(ii-1,jk,nibm,nitm)775 END DO776 END IF777 778 END IF ! End of array swap779 780 ! 2 - Calculation of radiation velocities781 ! ---------------------------------------782 783 IF( kt >= nit000 +3 .OR. ln_rstart ) THEN784 785 ! 2.1 Calculate the normal velocity based on phase velocity u_cynbnd786 ! -------------------------------------------------------------------787 !788 ! ji-row789 ! |790 ! nib -///u////// jpjnob + 1791 ! /////|//////792 ! nib -----f----- jpjnob793 ! |794 ! nibm-- u ---- jpjnob795 ! |796 ! nibm -----f----- jpjnob-1797 ! |798 ! nibm2-- u ---- jpjnob-1799 ! |800 ! nibm2 -----f----- jpjnob-2801 ! |802 ! ... radiative condition803 ! ... jpjnob+1,(jpindp1, jpinfm1)804 DO jj = fs_njn0+1, fs_njn1+1 ! Vector opt.805 DO jk = 1, jpkm1806 DO ji = 2, jpim1807 ! ... 2* j-gradient of u (f-point i=nibm, time mean)808 z2dx = ( unbnd(ji,jk,nibm ,nit) + unbnd(ji,jk,nibm ,nitm2) &809 - 2.*unbnd(ji,jk,nibm2,nitm)) / e2f(ji,jj-2)810 ! ... 2* i-gradient of u (u-point i=nibm, time nitm)811 z2dy = ( unbnd(ji+1,jk,nibm,nitm) - unbnd(ji-1,jk,nibm,nitm) ) / e1u(ji,jj-1)812 ! ... square of the norm of grad(v)813 z4nor2 = z2dx * z2dx + z2dy * z2dy814 ! ... minus time derivative (leap-frog) at nibm, without / 2 dt815 zdt = unbnd(ji,jk,nibm,nitm2) - unbnd(ji,jk,nibm,nit)816 ! ... i-phase speed ratio (bounded by 1) and save the unbounded phase817 ! velocity ratio no divided by e1f for the tracer radiation818 IF( z4nor2 == 0.) THEN819 z4nor2=.000001820 END IF821 z05cx = zdt * z2dx / z4nor2822 u_cynbnd(ji,jk) = z05cx *unmsk(ji,jk)823 END DO824 END DO825 END DO826 IF( lk_mpp ) CALL mppobc(u_cynbnd,jpind,jpinf,jpjnob+1,jpk,1,jpi, numout )827 828 ! ... extremeties njn0,njn1829 ii = jpind + 1 - nimpp830 IF( ii >= 2 .AND. ii < jpim1 ) THEN831 DO jk = 1, jpkm1832 u_cynbnd(ii,jk) = u_cynbnd(ii+1,jk)833 END DO834 END IF835 ii = jpinf + 1 - nimpp836 IF( ii >= 2 .AND. ii < jpim1 ) THEN837 DO jk = 1, jpkm1838 u_cynbnd(ii,jk) = u_cynbnd(ii-1,jk)839 END DO840 END IF841 842 ! 2.2 Calculate the normal velocity based on phase velocity v_cynbnd843 ! ------------------------------------------------------------------844 !845 ! ji-row ji-row846 ! |847 ! /////|/////////////////848 ! nib -----f----v----f---- jpjnob849 ! | |850 ! nib - u -- T -- u ---- jpjnob851 ! | |852 ! nibm -----f----v----f---- jpjnob-1853 ! | |854 ! nibm -- u -- T -- u --- jpjnob-1855 ! | |856 ! nibm2 -----f----v----f---- jpjnob-2857 ! | |858 ! ... Free surface formulation:859 ! ... radiative conditions on the total part + relaxation toward climatology860 ! ... jpjnob,(jpindp1, jpinfm1)861 DO jj = fs_njn0, fs_njn1 ! Vector opt.862 DO jk = 1, jpkm1863 DO ji = 2, jpim1864 ! ... 2* gradj(v) (T-point i=nibm, time mean)865 ii = ji -1 + nimpp866 z2dx = ( vnbnd(ji,jk,nibm ,nit) + vnbnd(ji,jk,nibm ,nitm2) &867 - 2.*vnbnd(ji,jk,nibm2,nitm)) / e2t(ji,jj-1)868 ! ... 2* gradi(v) (v-point i=nibm, time nitm)869 z2dy = ( vnbnd(ji+1,jk,nibm,nitm) - vnbnd(ji-1,jk,nibm,nitm) ) / e1v(ji,jj-1)870 ! ... square of the norm of grad(u)871 z4nor2 = z2dx * z2dx + z2dy * z2dy872 ! ... minus time derivative (leap-frog) at nibm, without / 2 dt873 zdt = vnbnd(ji,jk,nibm,nitm2) - vnbnd(ji,jk,nibm,nit)874 ! ... j-phase speed ratio (bounded by 1)875 IF( z4nor2 == 0. ) THEN876 z4nor2=.00001877 END IF878 z05cx = zdt * z2dx / z4nor2879 v_cynbnd(ji,jk)=z05cx *vnmsk(ji,jk)880 END DO881 END DO882 END DO883 884 END IF885 886 END SUBROUTINE obc_rad_north887 888 889 SUBROUTINE obc_rad_south ( kt )890 !!------------------------------------------------------------------------------891 !! *** SUBROUTINE obc_rad_south ***892 !!893 !! ** Purpose :894 !! Perform swap of arrays to calculate radiative phase speeds at the open895 !! south boundary and calculate those phase speeds if this OBC is not fixed.896 !! In case of fixed OBC, this subrountine is not called.897 !!898 !! History :899 !! ! 95-03 (J.-M. Molines) Original from SPEM900 !! ! 97-07 (G. Madec, J.-M. Molines) additions901 !! ! 97-12 (M. Imbard) Mpp adaptation902 !! ! 00-06 (J.-M. Molines)903 !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90904 !!------------------------------------------------------------------------------905 !! * Arguments906 INTEGER, INTENT( in ) :: kt907 908 !! * Local declarations909 INTEGER :: ii910 REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy911 REAL(wp) :: zvcb, zvcbm, zvcbm2912 !!------------------------------------------------------------------------------913 914 ! 1. Swap arrays before calculating radiative velocities915 ! ------------------------------------------------------916 917 ! 1.1 zonal velocity918 ! --------------------919 920 IF( kt > nit000 .OR. ln_rstart ) THEN921 922 ! ... advance in time (time filter, array swap)923 DO jk = 1, jpkm1924 DO ji = 1, jpi925 ! ... fields nitm2 <== nitm926 usbnd(ji,jk,nib ,nitm2) = usbnd(ji,jk,nib ,nitm)*usmsk(ji,jk)927 usbnd(ji,jk,nibm ,nitm2) = usbnd(ji,jk,nibm ,nitm)*usmsk(ji,jk)928 usbnd(ji,jk,nibm2,nitm2) = usbnd(ji,jk,nibm2,nitm)*usmsk(ji,jk)929 END DO930 END DO931 932 DO jj = fs_njs0, fs_njs1 ! Vector opt.933 DO jk = 1, jpkm1934 DO ji = 1, jpi935 usbnd(ji,jk,nib ,nitm) = usbnd(ji,jk,nib, nit)*usmsk(ji,jk)936 usbnd(ji,jk,nibm ,nitm) = usbnd(ji,jk,nibm ,nit)*usmsk(ji,jk)937 usbnd(ji,jk,nibm2,nitm) = usbnd(ji,jk,nibm2,nit)*usmsk(ji,jk)938 ! ... fields nit <== now (kt+1)939 usbnd(ji,jk,nib ,nit) = un(ji,jj ,jk)*usmsk(ji,jk)940 usbnd(ji,jk,nibm ,nit) = un(ji,jj+1,jk)*usmsk(ji,jk)941 usbnd(ji,jk,nibm2,nit) = un(ji,jj+2,jk)*usmsk(ji,jk)942 END DO943 END DO944 END DO945 IF( lk_mpp ) CALL mppobc(usbnd,jpisd,jpisf,jpjsob,jpk*3*3,1,jpi, numout )946 947 ! ... extremeties njs0,njs1948 ii = jpisd + 1 - nimpp949 IF( ii >= 2 .AND. ii < jpim1 ) THEN950 DO jk = 1, jpkm1951 usbnd(ii,jk,nibm,nitm) = usbnd(ii+1,jk,nibm,nitm)952 END DO953 END IF954 ii = jpisf + 1 - nimpp955 IF( ii >= 2 .AND. ii < jpim1 ) THEN956 DO jk = 1, jpkm1957 usbnd(ii,jk,nibm,nitm) = usbnd(ii-1,jk,nibm,nitm)958 END DO959 END IF960 961 ! 1.2 normal velocity962 ! -------------------963 964 !.. advance in time (time filter, array swap)965 DO jk = 1, jpkm1966 DO ji = 1, jpi967 ! ... fields nitm2 <== nitm968 vsbnd(ji,jk,nib ,nitm2) = vsbnd(ji,jk,nib ,nitm)*vsmsk(ji,jk)969 vsbnd(ji,jk,nibm ,nitm2) = vsbnd(ji,jk,nibm ,nitm)*vsmsk(ji,jk)970 END DO971 END DO972 973 DO jj = fs_njs0, fs_njs1 ! Vector opt.974 DO jk = 1, jpkm1975 DO ji = 1, jpi976 vsbnd(ji,jk,nib ,nitm) = vsbnd(ji,jk,nib, nit)*vsmsk(ji,jk)977 vsbnd(ji,jk,nibm ,nitm) = vsbnd(ji,jk,nibm ,nit)*vsmsk(ji,jk)978 vsbnd(ji,jk,nibm2,nitm) = vsbnd(ji,jk,nibm2,nit)*vsmsk(ji,jk)979 ! ... total or baroclinic velocity at b, bm and bm2980 zvcb = vn (ji,jj,jk)981 zvcbm = vn (ji,jj+1,jk)982 zvcbm2 = vn (ji,jj+2,jk)983 ! ... fields nit <== now (kt+1)984 vsbnd(ji,jk,nib ,nit) = zvcb *vsmsk(ji,jk)985 vsbnd(ji,jk,nibm ,nit) = zvcbm *vsmsk(ji,jk)986 vsbnd(ji,jk,nibm2,nit) = zvcbm2 *vsmsk(ji,jk)987 END DO988 END DO989 END DO990 IF( lk_mpp ) CALL mppobc(vsbnd,jpisd,jpisf,jpjsob,jpk*3*3,1,jpi, numout )991 992 ! ... extremeties njs0,njs1993 ii = jpisd + 1 - nimpp994 IF( ii >= 2 .AND. ii < jpim1 ) THEN995 DO jk = 1, jpkm1996 vsbnd(ii,jk,nibm,nitm) = vsbnd(ii+1,jk,nibm,nitm)997 END DO998 END IF999 ii = jpisf + 1 - nimpp1000 IF( ii >= 2 .AND. ii < jpim1 ) THEN1001 DO jk = 1, jpkm11002 vsbnd(ii,jk,nibm,nitm) = vsbnd(ii-1,jk,nibm,nitm)1003 END DO1004 END IF1005 1006 ! 1.3 Temperature and salinity1007 ! ----------------------------1008 1009 ! ... advance in time (time filter, array swap)1010 DO jk = 1, jpkm11011 DO ji = 1, jpi1012 ! ... fields nitm <== nit plus time filter at the boundary1013 tsbnd(ji,jk,nib,nitm) = tsbnd(ji,jk,nib,nit)*tsmsk(ji,jk)1014 ssbnd(ji,jk,nib,nitm) = ssbnd(ji,jk,nib,nit)*tsmsk(ji,jk)1015 END DO1016 END DO1017 1018 DO jj = fs_njs0, fs_njs1 ! Vector opt.1019 DO jk = 1, jpkm11020 DO ji = 1, jpi1021 tsbnd(ji,jk,nibm ,nitm) = tsbnd(ji,jk,nibm ,nit)*tsmsk(ji,jk)1022 ssbnd(ji,jk,nibm ,nitm) = ssbnd(ji,jk,nibm ,nit)*tsmsk(ji,jk)1023 ! ... fields nit <== now (kt+1)1024 tsbnd(ji,jk,nib ,nit) = tn(ji,jj ,jk)*tsmsk(ji,jk)1025 tsbnd(ji,jk,nibm ,nit) = tn(ji,jj+1 ,jk)*tsmsk(ji,jk)1026 ssbnd(ji,jk,nib ,nit) = sn(ji,jj ,jk)*tsmsk(ji,jk)1027 ssbnd(ji,jk,nibm ,nit) = sn(ji,jj+1 ,jk)*tsmsk(ji,jk)1028 END DO1029 END DO1030 END DO1031 IF( lk_mpp ) CALL mppobc(tsbnd,jpisd,jpisf,jpjsob,jpk*2*2,1,jpi, numout )1032 IF( lk_mpp ) CALL mppobc(ssbnd,jpisd,jpisf,jpjsob,jpk*2*2,1,jpi, numout )1033 1034 ! ... extremeties njs0,njs11035 ii = jpisd + 1 - nimpp1036 IF( ii >= 2 .AND. ii < jpim1 ) THEN1037 DO jk = 1, jpkm11038 tsbnd(ii,jk,nibm,nitm) = tsbnd(ii+1,jk,nibm,nitm)1039 ssbnd(ii,jk,nibm,nitm) = ssbnd(ii+1,jk,nibm,nitm)1040 END DO1041 END IF1042 ii = jpisf + 1 - nimpp1043 IF( ii >= 2 .AND. ii < jpim1 ) THEN1044 DO jk = 1, jpkm11045 tsbnd(ii,jk,nibm,nitm) = tsbnd(ii-1,jk,nibm,nitm)1046 ssbnd(ii,jk,nibm,nitm) = ssbnd(ii-1,jk,nibm,nitm)1047 END DO1048 END IF1049 1050 END IF ! End of array swap1051 1052 ! 2 - Calculation of radiation velocities1053 ! ---------------------------------------1054 1055 IF( kt >= nit000 +3 .OR. ln_rstart ) THEN1056 1057 ! 2.1 Calculate the normal velocity based on phase velocity u_cysbnd1058 ! -------------------------------------------------------------------1059 !1060 ! ji-row1061 ! |1062 ! nibm2 -----f----- jpjsob +21063 ! |1064 ! nibm2 -- u ----- jpjsob +21065 ! |1066 ! nibm -----f----- jpjsob +11067 ! |1068 ! nibm -- u ----- jpjsob +11069 ! |1070 ! nib -----f----- jpjsob1071 ! /////|//////1072 ! nib ////u///// jpjsob1073 !1074 ! ... radiative condition plus Raymond-Kuo1075 ! ... jpjsob,(jpisdp1, jpisfm1)1076 DO jj = fs_njs0, fs_njs1 ! Vector opt.1077 DO jk = 1, jpkm11078 DO ji = 2, jpim11079 ! ... 2* j-gradient of u (f-point i=nibm, time mean)1080 z2dx = (- usbnd(ji,jk,nibm ,nit) - usbnd(ji,jk,nibm ,nitm2) &1081 + 2.*usbnd(ji,jk,nibm2,nitm) ) / e2f(ji,jj+1)1082 ! ... 2* i-gradient of u (u-point i=nibm, time nitm)1083 z2dy = ( usbnd(ji+1,jk,nibm,nitm) - usbnd(ji-1,jk,nibm,nitm) ) / e1u(ji, jj+1)1084 ! ... square of the norm of grad(v)1085 z4nor2 = z2dx * z2dx + z2dy * z2dy1086 IF( z4nor2 == 0.) THEN1087 z4nor2 = 0.0000011088 END IF1089 ! ... minus time derivative (leap-frog) at nibm, without / 2 dt1090 zdt = usbnd(ji,jk,nibm,nitm2) - usbnd(ji,jk,nibm,nit)1091 ! ... i-phase speed ratio (bounded by -1) and save the unbounded phase1092 ! velocity ratio no divided by e1f for the tracer radiation1093 z05cx = zdt * z2dx / z4nor21094 u_cysbnd(ji,jk) = z05cx*usmsk(ji,jk)1095 END DO1096 END DO1097 END DO1098 IF( lk_mpp ) CALL mppobc(u_cysbnd,jpisd,jpisf,jpjsob,jpk,1,jpi, numout )1099 1100 ! ... extremeties njs0,njs11101 ii = jpisd + 1 - nimpp1102 IF( ii >= 2 .AND. ii < jpim1 ) THEN1103 DO jk = 1, jpkm11104 u_cysbnd(ii,jk) = u_cysbnd(ii+1,jk)1105 END DO1106 END IF1107 ii = jpisf + 1 - nimpp1108 IF( ii >= 2 .AND. ii < jpim1 ) THEN1109 DO jk = 1, jpkm11110 u_cysbnd(ii,jk) = u_cysbnd(ii-1,jk)1111 END DO1112 END IF1113 1114 ! 2.2 Calculate the normal velocity based on phase velocity v_cysbnd1115 ! -------------------------------------------------------------------1116 !1117 ! ji-row ji-row1118 ! | |1119 ! nibm2 -----f----v----f---- jpjsob+21120 ! | |1121 ! nibm - u -- T -- u ---- jpjsob+21122 ! | |1123 ! nibm -----f----v----f---- jpjsob+11124 ! | |1125 ! nib -- u -- T -- u --- jpjsob+11126 ! | |1127 ! nib -----f----v----f---- jpjsob1128 ! /////////////////////1129 !1130 ! ... Free surface formulation:1131 ! ... radiative conditions on the total part + relaxation toward climatology1132 ! ... jpjsob,(jpisdp1,jpisfm1)1133 DO jj = fs_njs0, fs_njs1 ! Vector opt.1134 DO jk = 1, jpkm11135 DO ji = 2, jpim11136 ! ... 2* gradj(v) (T-point i=nibm, time mean)1137 z2dx = ( - vsbnd(ji,jk,nibm ,nit) - vsbnd(ji,jk,nibm ,nitm2) &1138 + 2.*vsbnd(ji,jk,nibm2,nitm) ) / e2t(ji,jj+1)1139 ! ... 2* gradi(v) (v-point i=nibm, time nitm)1140 z2dy = ( vsbnd(ji+1,jk,nibm,nitm) - vsbnd(ji-1,jk,nibm,nitm) ) / e1v(ji,jj+1)1141 ! ... square of the norm of grad(u)1142 z4nor2 = z2dx * z2dx + z2dy * z2dy1143 IF( z4nor2 == 0.) THEN1144 z4nor2 = 0.0000011145 END IF1146 ! ... minus time derivative (leap-frog) at nibm, without / 2 dt1147 zdt = vsbnd(ji,jk,nibm,nitm2) - vsbnd(ji,jk,nibm,nit)1148 ! ... j-phase speed ratio (bounded by -1)1149 z05cx = zdt * z2dx / z4nor21150 v_cysbnd(ji,jk)=z05cx*vsmsk(ji,jk)1151 END DO1152 END DO1153 END DO1154 1155 ENDIF1156 1157 END SUBROUTINE obc_rad_south1158 1159 #else7 !!$ !!--------------------------------------------------------------------------------- 8 !!$ !! obc_rad : call the subroutine for each open boundary 9 !!$ !! obc_rad_east : compute the east phase velocities 10 !!$ !! obc_rad_west : compute the west phase velocities 11 !!$ !! obc_rad_north : compute the north phase velocities 12 !!$ !! obc_rad_south : compute the south phase velocities 13 !!$ !!--------------------------------------------------------------------------------- 14 !!$ USE oce ! ocean dynamics and tracers variables 15 !!$ USE dom_oce ! ocean space and time domain variables 16 !!$ USE lbclnk ! ocean lateral boundary conditions (or mpp link) 17 !!$ USE phycst ! physical constants 18 !!$ USE obc_oce ! ocean open boundary conditions 19 !!$ USE lib_mpp ! for mppobc 20 !!$ USE in_out_manager ! I/O units 21 !!$ 22 !!$ IMPLICIT NONE 23 !!$ PRIVATE 24 !!$ 25 !!$ PUBLIC obc_rad ! routine called by step.F90 26 !!$ 27 !!$ INTEGER :: ji, jj, jk ! dummy loop indices 28 !!$ 29 !!$ INTEGER :: & ! ... boundary space indices 30 !!$ nib = 1, & ! nib = boundary point 31 !!$ nibm = 2, & ! nibm = 1st interior point 32 !!$ nibm2 = 3, & ! nibm2 = 2nd interior point 33 !!$ ! ... boundary time indices 34 !!$ nit = 1, & ! nit = now 35 !!$ nitm = 2, & ! nitm = before 36 !!$ nitm2 = 3 ! nitm2 = before-before 37 !!$ 38 !!$ !! * Substitutions 39 !!$# include "obc_vectopt_loop_substitute.h90" 40 !!$ !!--------------------------------------------------------------------------------- 41 !!$ !! NEMO/OPA 3.3 , NEMO Consortium (2010) 42 !!$ !! $Id$ 43 !!$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 44 !!$ !!--------------------------------------------------------------------------------- 45 !!$ 46 !!$CONTAINS 47 !!$ 48 !!$ SUBROUTINE obc_rad ( kt ) 49 !!$ !!------------------------------------------------------------------------------ 50 !!$ !! SUBROUTINE obc_rad 51 !!$ !! ******************** 52 !!$ !! ** Purpose : 53 !!$ !! Perform swap of arrays to calculate radiative phase speeds at the open 54 !!$ !! boundaries and calculate those phase speeds if the open boundaries are 55 !!$ !! not fixed. In case of fixed open boundaries does nothing. 56 !!$ !! 57 !!$ !! The logical variable lp_obc_east, and/or lp_obc_west, and/or lp_obc_north, 58 !!$ !! and/or lp_obc_south allow the user to determine which boundary is an 59 !!$ !! open one (must be done in the param_obc.h90 file). 60 !!$ !! 61 !!$ !! ** Reference : 62 !!$ !! Marchesiello P., 1995, these de l'universite J. Fourier, Grenoble, France. 63 !!$ !! 64 !!$ !! History : 65 !!$ !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90 from the 66 !!$ !! J. Molines and G. Madec version 67 !!$ !!------------------------------------------------------------------------------ 68 !!$ INTEGER, INTENT( in ) :: kt 69 !!$ !!---------------------------------------------------------------------- 70 !!$ 71 !!$ IF( lp_obc_east .AND. .NOT.lfbceast ) CALL obc_rad_east ( kt ) ! East open boundary 72 !!$ 73 !!$ IF( lp_obc_west .AND. .NOT.lfbcwest ) CALL obc_rad_west ( kt ) ! West open boundary 74 !!$ 75 !!$ IF( lp_obc_north .AND. .NOT.lfbcnorth ) CALL obc_rad_north( kt ) ! North open boundary 76 !!$ 77 !!$ IF( lp_obc_south .AND. .NOT.lfbcsouth ) CALL obc_rad_south( kt ) ! South open boundary 78 !!$ 79 !!$ END SUBROUTINE obc_rad 80 !!$ 81 !!$ 82 !!$ SUBROUTINE obc_rad_east ( kt ) 83 !!$ !!------------------------------------------------------------------------------ 84 !!$ !! *** SUBROUTINE obc_rad_east *** 85 !!$ !! 86 !!$ !! ** Purpose : 87 !!$ !! Perform swap of arrays to calculate radiative phase speeds at the open 88 !!$ !! east boundary and calculate those phase speeds if this OBC is not fixed. 89 !!$ !! In case of fixed OBC, this subrountine is not called. 90 !!$ !! 91 !!$ !! History : 92 !!$ !! ! 95-03 (J.-M. Molines) Original from SPEM 93 !!$ !! ! 97-07 (G. Madec, J.-M. Molines) additions 94 !!$ !! ! 97-12 (M. Imbard) Mpp adaptation 95 !!$ !! ! 00-06 (J.-M. Molines) 96 !!$ !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90 97 !!$ !!------------------------------------------------------------------------------ 98 !!$ !! * Arguments 99 !!$ INTEGER, INTENT( in ) :: kt 100 !!$ 101 !!$ !! * Local declarations 102 !!$ INTEGER :: ij 103 !!$ REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy 104 !!$ REAL(wp) :: zucb, zucbm, zucbm2 105 !!$ !!------------------------------------------------------------------------------ 106 !!$ 107 !!$ ! 1. Swap arrays before calculating radiative velocities 108 !!$ ! ------------------------------------------------------ 109 !!$ 110 !!$ ! 1.1 zonal velocity 111 !!$ ! ------------------- 112 !!$ 113 !!$ IF( kt > nit000 .OR. ln_rstart ) THEN 114 !!$ 115 !!$ ! ... advance in time (time filter, array swap) 116 !!$ DO jk = 1, jpkm1 117 !!$ DO jj = 1, jpj 118 !!$ uebnd(jj,jk,nib ,nitm2) = uebnd(jj,jk,nib ,nitm)*uemsk(jj,jk) 119 !!$ uebnd(jj,jk,nibm ,nitm2) = uebnd(jj,jk,nibm ,nitm)*uemsk(jj,jk) 120 !!$ uebnd(jj,jk,nibm2,nitm2) = uebnd(jj,jk,nibm2,nitm)*uemsk(jj,jk) 121 !!$ END DO 122 !!$ END DO 123 !!$ ! ... fields nitm <== nit plus time filter at the boundary 124 !!$ DO ji = fs_nie0, fs_nie1 ! Vector opt. 125 !!$ DO jk = 1, jpkm1 126 !!$ DO jj = 1, jpj 127 !!$ uebnd(jj,jk,nib ,nitm) = uebnd(jj,jk,nib, nit)*uemsk(jj,jk) 128 !!$ uebnd(jj,jk,nibm ,nitm) = uebnd(jj,jk,nibm ,nit)*uemsk(jj,jk) 129 !!$ uebnd(jj,jk,nibm2,nitm) = uebnd(jj,jk,nibm2,nit)*uemsk(jj,jk) 130 !!$ ! ... fields nit <== now (kt+1) 131 !!$ ! ... Total or baroclinic velocity at b, bm and bm2 132 !!$ zucb = un(ji,jj,jk) 133 !!$ zucbm = un(ji-1,jj,jk) 134 !!$ zucbm2 = un(ji-2,jj,jk) 135 !!$ uebnd(jj,jk,nib ,nit) = zucb *uemsk(jj,jk) 136 !!$ uebnd(jj,jk,nibm ,nit) = zucbm *uemsk(jj,jk) 137 !!$ uebnd(jj,jk,nibm2,nit) = zucbm2 *uemsk(jj,jk) 138 !!$ END DO 139 !!$ END DO 140 !!$ END DO 141 !!$ IF( lk_mpp ) CALL mppobc(uebnd,jpjed,jpjef,jpieob,jpk*3*3,2,jpj, numout ) 142 !!$ 143 !!$ ! ... extremeties nie0, nie1 144 !!$ ij = jpjed +1 - njmpp 145 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 146 !!$ DO jk = 1,jpkm1 147 !!$ uebnd(ij,jk,nibm,nitm) = uebnd(ij+1 ,jk,nibm,nitm) 148 !!$ END DO 149 !!$ END IF 150 !!$ ij = jpjef +1 - njmpp 151 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 152 !!$ DO jk = 1,jpkm1 153 !!$ uebnd(ij,jk,nibm,nitm) = uebnd(ij-1 ,jk,nibm,nitm) 154 !!$ END DO 155 !!$ END IF 156 !!$ 157 !!$ ! 1.2 tangential velocity 158 !!$ ! ----------------------- 159 !!$ 160 !!$ ! ... advance in time (time filter, array swap) 161 !!$ DO jk = 1, jpkm1 162 !!$ DO jj = 1, jpj 163 !!$ ! ... fields nitm2 <== nitm 164 !!$ vebnd(jj,jk,nib ,nitm2) = vebnd(jj,jk,nib ,nitm)*vemsk(jj,jk) 165 !!$ vebnd(jj,jk,nibm ,nitm2) = vebnd(jj,jk,nibm ,nitm)*vemsk(jj,jk) 166 !!$ vebnd(jj,jk,nibm2,nitm2) = vebnd(jj,jk,nibm2,nitm)*vemsk(jj,jk) 167 !!$ END DO 168 !!$ END DO 169 !!$ 170 !!$ DO ji = fs_nie0+1, fs_nie1+1 ! Vector opt. 171 !!$ DO jk = 1, jpkm1 172 !!$ DO jj = 1, jpj 173 !!$ vebnd(jj,jk,nib ,nitm) = vebnd(jj,jk,nib, nit)*vemsk(jj,jk) 174 !!$ vebnd(jj,jk,nibm ,nitm) = vebnd(jj,jk,nibm ,nit)*vemsk(jj,jk) 175 !!$ vebnd(jj,jk,nibm2,nitm) = vebnd(jj,jk,nibm2,nit)*vemsk(jj,jk) 176 !!$ ! ... fields nit <== now (kt+1) 177 !!$ vebnd(jj,jk,nib ,nit) = vn(ji ,jj,jk)*vemsk(jj,jk) 178 !!$ vebnd(jj,jk,nibm ,nit) = vn(ji-1,jj,jk)*vemsk(jj,jk) 179 !!$ vebnd(jj,jk,nibm2,nit) = vn(ji-2,jj,jk)*vemsk(jj,jk) 180 !!$ END DO 181 !!$ END DO 182 !!$ END DO 183 !!$ IF( lk_mpp ) CALL mppobc(vebnd,jpjed,jpjef,jpieob+1,jpk*3*3,2,jpj, numout ) 184 !!$ 185 !!$ !... extremeties nie0, nie1 186 !!$ ij = jpjed +1 - njmpp 187 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 188 !!$ DO jk = 1,jpkm1 189 !!$ vebnd(ij,jk,nibm,nitm) = vebnd(ij+1 ,jk,nibm,nitm) 190 !!$ END DO 191 !!$ END IF 192 !!$ ij = jpjef +1 - njmpp 193 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 194 !!$ DO jk = 1,jpkm1 195 !!$ vebnd(ij,jk,nibm,nitm) = vebnd(ij-1 ,jk,nibm,nitm) 196 !!$ END DO 197 !!$ END IF 198 !!$ 199 !!$ ! 1.3 Temperature and salinity 200 !!$ ! ---------------------------- 201 !!$ 202 !!$ ! ... advance in time (time filter, array swap) 203 !!$ DO jk = 1, jpkm1 204 !!$ DO jj = 1, jpj 205 !!$ ! ... fields nitm <== nit plus time filter at the boundary 206 !!$ tebnd(jj,jk,nib,nitm) = tebnd(jj,jk,nib,nit)*temsk(jj,jk) 207 !!$ sebnd(jj,jk,nib,nitm) = sebnd(jj,jk,nib,nit)*temsk(jj,jk) 208 !!$ END DO 209 !!$ END DO 210 !!$ 211 !!$ DO ji = fs_nie0+1, fs_nie1+1 ! Vector opt. 212 !!$ DO jk = 1, jpkm1 213 !!$ DO jj = 1, jpj 214 !!$ tebnd(jj,jk,nibm,nitm) = tebnd(jj,jk,nibm,nit)*temsk(jj,jk) 215 !!$ sebnd(jj,jk,nibm,nitm) = sebnd(jj,jk,nibm,nit)*temsk(jj,jk) 216 !!$ ! ... fields nit <== now (kt+1) 217 !!$ tebnd(jj,jk,nib ,nit) = tn(ji ,jj,jk)*temsk(jj,jk) 218 !!$ tebnd(jj,jk,nibm ,nit) = tn(ji-1,jj,jk)*temsk(jj,jk) 219 !!$ sebnd(jj,jk,nib ,nit) = sn(ji ,jj,jk)*temsk(jj,jk) 220 !!$ sebnd(jj,jk,nibm ,nit) = sn(ji-1,jj,jk)*temsk(jj,jk) 221 !!$ END DO 222 !!$ END DO 223 !!$ END DO 224 !!$ IF( lk_mpp ) CALL mppobc(tebnd,jpjed,jpjef,jpieob+1,jpk*2*2,2,jpj, numout ) 225 !!$ IF( lk_mpp ) CALL mppobc(sebnd,jpjed,jpjef,jpieob+1,jpk*2*2,2,jpj, numout ) 226 !!$ 227 !!$ ! ... extremeties nie0, nie1 228 !!$ ij = jpjed +1 - njmpp 229 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 230 !!$ DO jk = 1,jpkm1 231 !!$ tebnd(ij,jk,nibm,nitm) = tebnd(ij+1 ,jk,nibm,nitm) 232 !!$ sebnd(ij,jk,nibm,nitm) = sebnd(ij+1 ,jk,nibm,nitm) 233 !!$ END DO 234 !!$ END IF 235 !!$ ij = jpjef +1 - njmpp 236 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 237 !!$ DO jk = 1,jpkm1 238 !!$ tebnd(ij,jk,nibm,nitm) = tebnd(ij-1 ,jk,nibm,nitm) 239 !!$ sebnd(ij,jk,nibm,nitm) = sebnd(ij-1 ,jk,nibm,nitm) 240 !!$ END DO 241 !!$ END IF 242 !!$ 243 !!$ END IF ! End of array swap 244 !!$ 245 !!$ ! 2 - Calculation of radiation velocities 246 !!$ ! --------------------------------------- 247 !!$ 248 !!$ IF( kt >= nit000 +3 .OR. ln_rstart ) THEN 249 !!$ 250 !!$ ! 2.1 Calculate the normal velocity U based on phase velocity u_cxebnd 251 !!$ ! --------------------------------------------------------------------- 252 !!$ ! 253 !!$ ! nibm2 nibm nib 254 !!$ ! | nibm | nib |/// 255 !!$ ! | | | | |/// 256 !!$ ! jj-line --f----v----f----v----f--- 257 !!$ ! | | | | |/// 258 !!$ ! | | |/// 259 !!$ ! jj-line u T u T u/// 260 !!$ ! | | |/// 261 !!$ ! | | | | |/// 262 !!$ ! jpieob-2 jpieob-1 jpieob 263 !!$ ! | | 264 !!$ ! jpieob-1 jpieob 265 !!$ ! 266 !!$ ! ... (jpjedp1, jpjefm1),jpieob 267 !!$ DO ji = fs_nie0, fs_nie1 ! Vector opt. 268 !!$ DO jk = 1, jpkm1 269 !!$ DO jj = 2, jpjm1 270 !!$ ! ... 2* gradi(u) (T-point i=nibm, time mean) 271 !!$ z2dx = ( uebnd(jj,jk,nibm ,nit) + uebnd(jj,jk,nibm ,nitm2) & 272 !!$ - 2.*uebnd(jj,jk,nibm2,nitm) ) / e1t(ji-1,jj) 273 !!$ ! ... 2* gradj(u) (u-point i=nibm, time nitm) 274 !!$ z2dy = ( uebnd(jj+1,jk,nibm,nitm) - uebnd(jj-1,jk,nibm,nitm) ) / e2u(ji-1,jj) 275 !!$ ! ... square of the norm of grad(u) 276 !!$ z4nor2 = z2dx * z2dx + z2dy * z2dy 277 !!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt 278 !!$ zdt = uebnd(jj,jk,nibm,nitm2) - uebnd(jj,jk,nibm,nit) 279 !!$ ! ... i-phase speed ratio (bounded by 1) 280 !!$ IF( z4nor2 == 0. ) THEN 281 !!$ z4nor2=.00001 282 !!$ END IF 283 !!$ z05cx = zdt * z2dx / z4nor2 284 !!$ u_cxebnd(jj,jk) = z05cx*uemsk(jj,jk) 285 !!$ END DO 286 !!$ END DO 287 !!$ END DO 288 !!$ 289 !!$ ! 2.2 Calculate the tangential velocity based on phase velocity v_cxebnd 290 !!$ ! ----------------------------------------------------------------------- 291 !!$ ! 292 !!$ ! nibm2 nibm nib 293 !!$ ! | nibm | nib///|/// 294 !!$ ! | | | |////|/// 295 !!$ ! jj-line --v----f----v----f----v--- 296 !!$ ! | | | |////|/// 297 !!$ ! | | | |////|/// 298 !!$ ! | jpieob-1| jpieob /|/// 299 !!$ ! | | | 300 !!$ ! jpieob-1 jpieob jpieob+1 301 !!$ ! 302 !!$ ! ... (jpjedp1, jpjefm1), jpieob+1 303 !!$ DO ji = fs_nie0+1, fs_nie1+1 ! Vector opt. 304 !!$ DO jk = 1, jpkm1 305 !!$ DO jj = 2, jpjm1 306 !!$ ! ... 2* i-gradient of v (f-point i=nibm, time mean) 307 !!$ z2dx = ( vebnd(jj,jk,nibm ,nit) + vebnd(jj,jk,nibm ,nitm2) & 308 !!$ - 2.*vebnd(jj,jk,nibm2,nitm) ) / e1f(ji-2,jj) 309 !!$ ! ... 2* j-gradient of v (v-point i=nibm, time nitm) 310 !!$ z2dy = ( vebnd(jj+1,jk,nibm,nitm) - vebnd(jj-1,jk,nibm,nitm) ) / e2v(ji-1,jj) 311 !!$ ! ... square of the norm of grad(v) 312 !!$ z4nor2 = z2dx * z2dx + z2dy * z2dy 313 !!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt 314 !!$ zdt = vebnd(jj,jk,nibm,nitm2) - vebnd(jj,jk,nibm,nit) 315 !!$ ! ... i-phase speed ratio (bounded by 1) and save the unbounded phase 316 !!$ ! velocity ratio no divided by e1f for the tracer radiation 317 !!$ IF( z4nor2 == 0. ) THEN 318 !!$ z4nor2=.000001 319 !!$ END IF 320 !!$ z05cx = zdt * z2dx / z4nor2 321 !!$ v_cxebnd(jj,jk) = z05cx*vemsk(jj,jk) 322 !!$ END DO 323 !!$ END DO 324 !!$ END DO 325 !!$ IF( lk_mpp ) CALL mppobc(v_cxebnd,jpjed,jpjef,jpieob+1,jpk,2,jpj, numout ) 326 !!$ 327 !!$ ! ... extremeties nie0, nie1 328 !!$ ij = jpjed +1 - njmpp 329 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 330 !!$ DO jk = 1,jpkm1 331 !!$ v_cxebnd(ij,jk) = v_cxebnd(ij+1 ,jk) 332 !!$ END DO 333 !!$ END IF 334 !!$ ij = jpjef +1 - njmpp 335 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 336 !!$ DO jk = 1,jpkm1 337 !!$ v_cxebnd(ij,jk) = v_cxebnd(ij-1 ,jk) 338 !!$ END DO 339 !!$ END IF 340 !!$ 341 !!$ END IF 342 !!$ 343 !!$ END SUBROUTINE obc_rad_east 344 !!$ 345 !!$ 346 !!$ SUBROUTINE obc_rad_west ( kt ) 347 !!$ !!------------------------------------------------------------------------------ 348 !!$ !! *** SUBROUTINE obc_rad_west *** 349 !!$ !! 350 !!$ !! ** Purpose : 351 !!$ !! Perform swap of arrays to calculate radiative phase speeds at the open 352 !!$ !! west boundary and calculate those phase speeds if this OBC is not fixed. 353 !!$ !! In case of fixed OBC, this subrountine is not called. 354 !!$ !! 355 !!$ !! History : 356 !!$ !! ! 95-03 (J.-M. Molines) Original from SPEM 357 !!$ !! ! 97-07 (G. Madec, J.-M. Molines) additions 358 !!$ !! ! 97-12 (M. Imbard) Mpp adaptation 359 !!$ !! ! 00-06 (J.-M. Molines) 360 !!$ !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90 361 !!$ !!------------------------------------------------------------------------------ 362 !!$ !! * Arguments 363 !!$ INTEGER, INTENT( in ) :: kt 364 !!$ 365 !!$ !! * Local declarations 366 !!$ INTEGER :: ij 367 !!$ REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy 368 !!$ REAL(wp) :: zucb, zucbm, zucbm2 369 !!$ !!------------------------------------------------------------------------------ 370 !!$ 371 !!$ ! 1. Swap arrays before calculating radiative velocities 372 !!$ ! ------------------------------------------------------ 373 !!$ 374 !!$ ! 1.1 zonal velocity 375 !!$ ! ------------------- 376 !!$ 377 !!$ IF( kt > nit000 .OR. ln_rstart ) THEN 378 !!$ 379 !!$ ! ... advance in time (time filter, array swap) 380 !!$ DO jk = 1, jpkm1 381 !!$ DO jj = 1, jpj 382 !!$ uwbnd(jj,jk,nib ,nitm2) = uwbnd(jj,jk,nib ,nitm)*uwmsk(jj,jk) 383 !!$ uwbnd(jj,jk,nibm ,nitm2) = uwbnd(jj,jk,nibm ,nitm)*uwmsk(jj,jk) 384 !!$ uwbnd(jj,jk,nibm2,nitm2) = uwbnd(jj,jk,nibm2,nitm)*uwmsk(jj,jk) 385 !!$ END DO 386 !!$ END DO 387 !!$ 388 !!$ ! ... fields nitm <== nit plus time filter at the boundary 389 !!$ DO ji = fs_niw0, fs_niw1 ! Vector opt. 390 !!$ DO jk = 1, jpkm1 391 !!$ DO jj = 1, jpj 392 !!$ uwbnd(jj,jk,nib ,nitm) = uwbnd(jj,jk,nib ,nit)*uwmsk(jj,jk) 393 !!$ uwbnd(jj,jk,nibm ,nitm) = uwbnd(jj,jk,nibm ,nit)*uwmsk(jj,jk) 394 !!$ uwbnd(jj,jk,nibm2,nitm) = uwbnd(jj,jk,nibm2,nit)*uwmsk(jj,jk) 395 !!$ ! ... total or baroclinic velocity at b, bm and bm2 396 !!$ zucb = un (ji,jj,jk) 397 !!$ zucbm = un (ji+1,jj,jk) 398 !!$ zucbm2 = un (ji+2,jj,jk) 399 !!$ 400 !!$ ! ... fields nit <== now (kt+1) 401 !!$ uwbnd(jj,jk,nib ,nit) = zucb *uwmsk(jj,jk) 402 !!$ uwbnd(jj,jk,nibm ,nit) = zucbm *uwmsk(jj,jk) 403 !!$ uwbnd(jj,jk,nibm2,nit) = zucbm2*uwmsk(jj,jk) 404 !!$ END DO 405 !!$ END DO 406 !!$ END DO 407 !!$ IF( lk_mpp ) CALL mppobc(uwbnd,jpjwd,jpjwf,jpiwob,jpk*3*3,2,jpj, numout ) 408 !!$ 409 !!$ ! ... extremeties niw0, niw1 410 !!$ ij = jpjwd +1 - njmpp 411 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 412 !!$ DO jk = 1,jpkm1 413 !!$ uwbnd(ij,jk,nibm,nitm) = uwbnd(ij+1 ,jk,nibm,nitm) 414 !!$ END DO 415 !!$ END IF 416 !!$ ij = jpjwf +1 - njmpp 417 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 418 !!$ DO jk = 1,jpkm1 419 !!$ uwbnd(ij,jk,nibm,nitm) = uwbnd(ij-1 ,jk,nibm,nitm) 420 !!$ END DO 421 !!$ END IF 422 !!$ 423 !!$ ! 1.2 tangential velocity 424 !!$ ! ----------------------- 425 !!$ 426 !!$ ! ... advance in time (time filter, array swap) 427 !!$ DO jk = 1, jpkm1 428 !!$ DO jj = 1, jpj 429 !!$ ! ... fields nitm2 <== nitm 430 !!$ vwbnd(jj,jk,nib ,nitm2) = vwbnd(jj,jk,nib ,nitm)*vwmsk(jj,jk) 431 !!$ vwbnd(jj,jk,nibm ,nitm2) = vwbnd(jj,jk,nibm ,nitm)*vwmsk(jj,jk) 432 !!$ vwbnd(jj,jk,nibm2,nitm2) = vwbnd(jj,jk,nibm2,nitm)*vwmsk(jj,jk) 433 !!$ END DO 434 !!$ END DO 435 !!$ 436 !!$ DO ji = fs_niw0, fs_niw1 ! Vector opt. 437 !!$ DO jk = 1, jpkm1 438 !!$ DO jj = 1, jpj 439 !!$ vwbnd(jj,jk,nib ,nitm) = vwbnd(jj,jk,nib, nit)*vwmsk(jj,jk) 440 !!$ vwbnd(jj,jk,nibm ,nitm) = vwbnd(jj,jk,nibm ,nit)*vwmsk(jj,jk) 441 !!$ vwbnd(jj,jk,nibm2,nitm) = vwbnd(jj,jk,nibm2,nit)*vwmsk(jj,jk) 442 !!$ ! ... fields nit <== now (kt+1) 443 !!$ vwbnd(jj,jk,nib ,nit) = vn(ji ,jj,jk)*vwmsk(jj,jk) 444 !!$ vwbnd(jj,jk,nibm ,nit) = vn(ji+1,jj,jk)*vwmsk(jj,jk) 445 !!$ vwbnd(jj,jk,nibm2,nit) = vn(ji+2,jj,jk)*vwmsk(jj,jk) 446 !!$ END DO 447 !!$ END DO 448 !!$ END DO 449 !!$ IF( lk_mpp ) CALL mppobc(vwbnd,jpjwd,jpjwf,jpiwob,jpk*3*3,2,jpj, numout ) 450 !!$ 451 !!$ ! ... extremeties niw0, niw1 452 !!$ ij = jpjwd +1 - njmpp 453 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 454 !!$ DO jk = 1,jpkm1 455 !!$ vwbnd(ij,jk,nibm,nitm) = vwbnd(ij+1 ,jk,nibm,nitm) 456 !!$ END DO 457 !!$ END IF 458 !!$ ij = jpjwf +1 - njmpp 459 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 460 !!$ DO jk = 1,jpkm1 461 !!$ vwbnd(ij,jk,nibm,nitm) = vwbnd(ij-1 ,jk,nibm,nitm) 462 !!$ END DO 463 !!$ END IF 464 !!$ 465 !!$ ! 1.3 Temperature and salinity 466 !!$ ! ---------------------------- 467 !!$ 468 !!$ ! ... advance in time (time filter, array swap) 469 !!$ DO jk = 1, jpkm1 470 !!$ DO jj = 1, jpj 471 !!$ ! ... fields nitm <== nit plus time filter at the boundary 472 !!$ twbnd(jj,jk,nib,nitm) = twbnd(jj,jk,nib,nit)*twmsk(jj,jk) 473 !!$ swbnd(jj,jk,nib,nitm) = swbnd(jj,jk,nib,nit)*twmsk(jj,jk) 474 !!$ END DO 475 !!$ END DO 476 !!$ 477 !!$ DO ji = fs_niw0, fs_niw1 ! Vector opt. 478 !!$ DO jk = 1, jpkm1 479 !!$ DO jj = 1, jpj 480 !!$ twbnd(jj,jk,nibm ,nitm) = twbnd(jj,jk,nibm ,nit)*twmsk(jj,jk) 481 !!$ swbnd(jj,jk,nibm ,nitm) = swbnd(jj,jk,nibm ,nit)*twmsk(jj,jk) 482 !!$ ! ... fields nit <== now (kt+1) 483 !!$ twbnd(jj,jk,nib ,nit) = tn(ji ,jj,jk)*twmsk(jj,jk) 484 !!$ twbnd(jj,jk,nibm ,nit) = tn(ji+1 ,jj,jk)*twmsk(jj,jk) 485 !!$ swbnd(jj,jk,nib ,nit) = sn(ji ,jj,jk)*twmsk(jj,jk) 486 !!$ swbnd(jj,jk,nibm ,nit) = sn(ji+1 ,jj,jk)*twmsk(jj,jk) 487 !!$ END DO 488 !!$ END DO 489 !!$ END DO 490 !!$ IF( lk_mpp ) CALL mppobc(twbnd,jpjwd,jpjwf,jpiwob,jpk*2*2,2,jpj, numout ) 491 !!$ IF( lk_mpp ) CALL mppobc(swbnd,jpjwd,jpjwf,jpiwob,jpk*2*2,2,jpj, numout ) 492 !!$ 493 !!$ ! ... extremeties niw0, niw1 494 !!$ ij = jpjwd +1 - njmpp 495 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 496 !!$ DO jk = 1,jpkm1 497 !!$ twbnd(ij,jk,nibm,nitm) = twbnd(ij+1 ,jk,nibm,nitm) 498 !!$ swbnd(ij,jk,nibm,nitm) = swbnd(ij+1 ,jk,nibm,nitm) 499 !!$ END DO 500 !!$ END IF 501 !!$ ij = jpjwf +1 - njmpp 502 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 503 !!$ DO jk = 1,jpkm1 504 !!$ twbnd(ij,jk,nibm,nitm) = twbnd(ij-1 ,jk,nibm,nitm) 505 !!$ swbnd(ij,jk,nibm,nitm) = swbnd(ij-1 ,jk,nibm,nitm) 506 !!$ END DO 507 !!$ END IF 508 !!$ 509 !!$ END IF ! End of array swap 510 !!$ 511 !!$ ! 2 - Calculation of radiation velocities 512 !!$ ! --------------------------------------- 513 !!$ 514 !!$ IF( kt >= nit000 +3 .OR. ln_rstart ) THEN 515 !!$ 516 !!$ ! 2.1 Calculate the normal velocity U based on phase velocity u_cxwbnd 517 !!$ ! ---------------------------------------------------------------------- 518 !!$ ! 519 !!$ ! nib nibm nibm2 520 !!$ ! ///| nib | nibm | 521 !!$ ! ///| | | | | 522 !!$ ! ---f----v----f----v----f-- jj-line 523 !!$ ! ///| | | | | 524 !!$ ! ///| | | 525 !!$ ! ///u T u T u jj-line 526 !!$ ! ///| | | 527 !!$ ! ///| | | | | 528 !!$ ! jpiwob jpiwob+1 jpiwob+2 529 !!$ ! | | 530 !!$ ! jpiwob+1 jpiwob+2 531 !!$ ! 532 !!$ ! ... If free surface formulation: 533 !!$ ! ... radiative conditions on the total part + relaxation toward climatology 534 !!$ ! ... (jpjwdp1, jpjwfm1), jpiwob 535 !!$ DO ji = fs_niw0, fs_niw1 ! Vector opt. 536 !!$ DO jk = 1, jpkm1 537 !!$ DO jj = 2, jpjm1 538 !!$ ! ... 2* gradi(u) (T-point i=nibm, time mean) 539 !!$ z2dx = ( - uwbnd(jj,jk,nibm ,nit) - uwbnd(jj,jk,nibm ,nitm2) & 540 !!$ + 2.*uwbnd(jj,jk,nibm2,nitm) ) / e1t(ji+2,jj) 541 !!$ ! ... 2* gradj(u) (u-point i=nibm, time nitm) 542 !!$ z2dy = ( uwbnd(jj+1,jk,nibm,nitm) - uwbnd(jj-1,jk,nibm,nitm) ) / e2u(ji+1,jj) 543 !!$ ! ... square of the norm of grad(u) 544 !!$ z4nor2 = z2dx * z2dx + z2dy * z2dy 545 !!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt 546 !!$ zdt = uwbnd(jj,jk,nibm,nitm2) - uwbnd(jj,jk,nibm,nit) 547 !!$ ! ... i-phase speed ratio (bounded by -1) 548 !!$ IF( z4nor2 == 0. ) THEN 549 !!$ z4nor2=0.00001 550 !!$ END IF 551 !!$ z05cx = zdt * z2dx / z4nor2 552 !!$ u_cxwbnd(jj,jk)=z05cx*uwmsk(jj,jk) 553 !!$ END DO 554 !!$ END DO 555 !!$ END DO 556 !!$ 557 !!$ ! 2.2 Calculate the tangential velocity based on phase velocity v_cxwbnd 558 !!$ ! ----------------------------------------------------------------------- 559 !!$ ! 560 !!$ ! nib nibm nibm2 561 !!$ ! ///|///nib | nibm | nibm2 562 !!$ ! ///|////| | | | | | 563 !!$ ! ---v----f----v----f----v----f----v-- jj-line 564 !!$ ! ///|////| | | | | | 565 !!$ ! ///|////| | | | | | 566 !!$ ! jpiwob jpiwob+1 jpiwob+2 567 !!$ ! | | | 568 !!$ ! jpiwob jpiwob+1 jpiwob+2 569 !!$ ! 570 !!$ ! ... radiative condition plus Raymond-Kuo 571 !!$ ! ... (jpjwdp1, jpjwfm1),jpiwob 572 !!$ DO ji = fs_niw0, fs_niw1 ! Vector opt. 573 !!$ DO jk = 1, jpkm1 574 !!$ DO jj = 2, jpjm1 575 !!$ ! ... 2* i-gradient of v (f-point i=nibm, time mean) 576 !!$ z2dx = ( - vwbnd(jj,jk,nibm ,nit) - vwbnd(jj,jk,nibm ,nitm2) & 577 !!$ + 2.*vwbnd(jj,jk,nibm2,nitm) ) / e1f(ji+1,jj) 578 !!$ ! ... 2* j-gradient of v (v-point i=nibm, time nitm) 579 !!$ z2dy = ( vwbnd(jj+1,jk,nibm,nitm) - vwbnd(jj-1,jk,nibm,nitm) ) / e2v(ji+1,jj) 580 !!$ ! ... square of the norm of grad(v) 581 !!$ z4nor2 = z2dx * z2dx + z2dy * z2dy 582 !!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt 583 !!$ zdt = vwbnd(jj,jk,nibm,nitm2) - vwbnd(jj,jk,nibm,nit) 584 !!$ ! ... i-phase speed ratio (bounded by -1) and save the unbounded phase 585 !!$ ! velocity ratio no divided by e1f for the tracer radiation 586 !!$ IF( z4nor2 == 0) THEN 587 !!$ z4nor2=0.000001 588 !!$ endif 589 !!$ z05cx = zdt * z2dx / z4nor2 590 !!$ v_cxwbnd(jj,jk) = z05cx*vwmsk(jj,jk) 591 !!$ END DO 592 !!$ END DO 593 !!$ END DO 594 !!$ IF( lk_mpp ) CALL mppobc(v_cxwbnd,jpjwd,jpjwf,jpiwob,jpk,2,jpj, numout ) 595 !!$ 596 !!$ ! ... extremeties niw0, niw1 597 !!$ ij = jpjwd +1 - njmpp 598 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 599 !!$ DO jk = 1,jpkm1 600 !!$ v_cxwbnd(ij,jk) = v_cxwbnd(ij+1 ,jk) 601 !!$ END DO 602 !!$ END IF 603 !!$ ij = jpjwf +1 - njmpp 604 !!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN 605 !!$ DO jk = 1,jpkm1 606 !!$ v_cxwbnd(ij,jk) = v_cxwbnd(ij-1 ,jk) 607 !!$ END DO 608 !!$ END IF 609 !!$ 610 !!$ END IF 611 !!$ 612 !!$ END SUBROUTINE obc_rad_west 613 !!$ 614 !!$ 615 !!$ SUBROUTINE obc_rad_north ( kt ) 616 !!$ !!------------------------------------------------------------------------------ 617 !!$ !! *** SUBROUTINE obc_rad_north *** 618 !!$ !! 619 !!$ !! ** Purpose : 620 !!$ !! Perform swap of arrays to calculate radiative phase speeds at the open 621 !!$ !! north boundary and calculate those phase speeds if this OBC is not fixed. 622 !!$ !! In case of fixed OBC, this subrountine is not called. 623 !!$ !! 624 !!$ !! History : 625 !!$ !! ! 95-03 (J.-M. Molines) Original from SPEM 626 !!$ !! ! 97-07 (G. Madec, J.-M. Molines) additions 627 !!$ !! ! 97-12 (M. Imbard) Mpp adaptation 628 !!$ !! ! 00-06 (J.-M. Molines) 629 !!$ !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90 630 !!$ !!------------------------------------------------------------------------------ 631 !!$ !! * Arguments 632 !!$ INTEGER, INTENT( in ) :: kt 633 !!$ 634 !!$ !! * Local declarations 635 !!$ INTEGER :: ii 636 !!$ REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy 637 !!$ REAL(wp) :: zvcb, zvcbm, zvcbm2 638 !!$ !!------------------------------------------------------------------------------ 639 !!$ 640 !!$ ! 1. Swap arrays before calculating radiative velocities 641 !!$ ! ------------------------------------------------------ 642 !!$ 643 !!$ ! 1.1 zonal velocity 644 !!$ ! ------------------- 645 !!$ 646 !!$ IF( kt > nit000 .OR. ln_rstart ) THEN 647 !!$ 648 !!$ ! ... advance in time (time filter, array swap) 649 !!$ DO jk = 1, jpkm1 650 !!$ DO ji = 1, jpi 651 !!$ ! ... fields nitm2 <== nitm 652 !!$ unbnd(ji,jk,nib ,nitm2) = unbnd(ji,jk,nib ,nitm)*unmsk(ji,jk) 653 !!$ unbnd(ji,jk,nibm ,nitm2) = unbnd(ji,jk,nibm ,nitm)*unmsk(ji,jk) 654 !!$ unbnd(ji,jk,nibm2,nitm2) = unbnd(ji,jk,nibm2,nitm)*unmsk(ji,jk) 655 !!$ END DO 656 !!$ END DO 657 !!$ 658 !!$ DO jj = fs_njn0+1, fs_njn1+1 ! Vector opt. 659 !!$ DO jk = 1, jpkm1 660 !!$ DO ji = 1, jpi 661 !!$ unbnd(ji,jk,nib ,nitm) = unbnd(ji,jk,nib, nit)*unmsk(ji,jk) 662 !!$ unbnd(ji,jk,nibm ,nitm) = unbnd(ji,jk,nibm ,nit)*unmsk(ji,jk) 663 !!$ unbnd(ji,jk,nibm2,nitm) = unbnd(ji,jk,nibm2,nit)*unmsk(ji,jk) 664 !!$ ! ... fields nit <== now (kt+1) 665 !!$ unbnd(ji,jk,nib ,nit) = un(ji,jj, jk)*unmsk(ji,jk) 666 !!$ unbnd(ji,jk,nibm ,nit) = un(ji,jj-1,jk)*unmsk(ji,jk) 667 !!$ unbnd(ji,jk,nibm2,nit) = un(ji,jj-2,jk)*unmsk(ji,jk) 668 !!$ END DO 669 !!$ END DO 670 !!$ END DO 671 !!$ IF( lk_mpp ) CALL mppobc(unbnd,jpind,jpinf,jpjnob+1,jpk*3*3,1,jpi, numout ) 672 !!$ 673 !!$ ! ... extremeties njn0,njn1 674 !!$ ii = jpind + 1 - nimpp 675 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 676 !!$ DO jk = 1, jpkm1 677 !!$ unbnd(ii,jk,nibm,nitm) = unbnd(ii+1,jk,nibm,nitm) 678 !!$ END DO 679 !!$ END IF 680 !!$ ii = jpinf + 1 - nimpp 681 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 682 !!$ DO jk = 1, jpkm1 683 !!$ unbnd(ii,jk,nibm,nitm) = unbnd(ii-1,jk,nibm,nitm) 684 !!$ END DO 685 !!$ END IF 686 !!$ 687 !!$ ! 1.2. normal velocity 688 !!$ ! -------------------- 689 !!$ 690 !!$ ! ... advance in time (time filter, array swap) 691 !!$ DO jk = 1, jpkm1 692 !!$ DO ji = 1, jpi 693 !!$ ! ... fields nitm2 <== nitm 694 !!$ vnbnd(ji,jk,nib ,nitm2) = vnbnd(ji,jk,nib ,nitm)*vnmsk(ji,jk) 695 !!$ vnbnd(ji,jk,nibm ,nitm2) = vnbnd(ji,jk,nibm ,nitm)*vnmsk(ji,jk) 696 !!$ vnbnd(ji,jk,nibm2,nitm2) = vnbnd(ji,jk,nibm2,nitm)*vnmsk(ji,jk) 697 !!$ END DO 698 !!$ END DO 699 !!$ 700 !!$ DO jj = fs_njn0, fs_njn1 ! Vector opt. 701 !!$ DO jk = 1, jpkm1 702 !!$ DO ji = 1, jpi 703 !!$ vnbnd(ji,jk,nib ,nitm) = vnbnd(ji,jk,nib, nit)*vnmsk(ji,jk) 704 !!$ vnbnd(ji,jk,nibm ,nitm) = vnbnd(ji,jk,nibm ,nit)*vnmsk(ji,jk) 705 !!$ vnbnd(ji,jk,nibm2,nitm) = vnbnd(ji,jk,nibm2,nit)*vnmsk(ji,jk) 706 !!$ ! ... fields nit <== now (kt+1) 707 !!$ ! ... total or baroclinic velocity at b, bm and bm2 708 !!$ zvcb = vn (ji,jj,jk) 709 !!$ zvcbm = vn (ji,jj-1,jk) 710 !!$ zvcbm2 = vn (ji,jj-2,jk) 711 !!$ ! ... fields nit <== now (kt+1) 712 !!$ vnbnd(ji,jk,nib ,nit) = zvcb *vnmsk(ji,jk) 713 !!$ vnbnd(ji,jk,nibm ,nit) = zvcbm *vnmsk(ji,jk) 714 !!$ vnbnd(ji,jk,nibm2,nit) = zvcbm2*vnmsk(ji,jk) 715 !!$ END DO 716 !!$ END DO 717 !!$ END DO 718 !!$ IF( lk_mpp ) CALL mppobc(vnbnd,jpind,jpinf,jpjnob,jpk*3*3,1,jpi, numout ) 719 !!$ 720 !!$ ! ... extremeties njn0,njn1 721 !!$ ii = jpind + 1 - nimpp 722 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 723 !!$ DO jk = 1, jpkm1 724 !!$ vnbnd(ii,jk,nibm,nitm) = vnbnd(ii+1,jk,nibm,nitm) 725 !!$ END DO 726 !!$ END IF 727 !!$ ii = jpinf + 1 - nimpp 728 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 729 !!$ DO jk = 1, jpkm1 730 !!$ vnbnd(ii,jk,nibm,nitm) = vnbnd(ii-1,jk,nibm,nitm) 731 !!$ END DO 732 !!$ END IF 733 !!$ 734 !!$ ! 1.3 Temperature and salinity 735 !!$ ! ---------------------------- 736 !!$ 737 !!$ ! ... advance in time (time filter, array swap) 738 !!$ DO jk = 1, jpkm1 739 !!$ DO ji = 1, jpi 740 !!$ ! ... fields nitm <== nit plus time filter at the boundary 741 !!$ tnbnd(ji,jk,nib ,nitm) = tnbnd(ji,jk,nib,nit)*tnmsk(ji,jk) 742 !!$ snbnd(ji,jk,nib ,nitm) = snbnd(ji,jk,nib,nit)*tnmsk(ji,jk) 743 !!$ END DO 744 !!$ END DO 745 !!$ 746 !!$ DO jj = fs_njn0+1, fs_njn1+1 ! Vector opt. 747 !!$ DO jk = 1, jpkm1 748 !!$ DO ji = 1, jpi 749 !!$ tnbnd(ji,jk,nibm ,nitm) = tnbnd(ji,jk,nibm ,nit)*tnmsk(ji,jk) 750 !!$ snbnd(ji,jk,nibm ,nitm) = snbnd(ji,jk,nibm ,nit)*tnmsk(ji,jk) 751 !!$ ! ... fields nit <== now (kt+1) 752 !!$ tnbnd(ji,jk,nib ,nit) = tn(ji,jj, jk)*tnmsk(ji,jk) 753 !!$ tnbnd(ji,jk,nibm ,nit) = tn(ji,jj-1,jk)*tnmsk(ji,jk) 754 !!$ snbnd(ji,jk,nib ,nit) = sn(ji,jj, jk)*tnmsk(ji,jk) 755 !!$ snbnd(ji,jk,nibm ,nit) = sn(ji,jj-1,jk)*tnmsk(ji,jk) 756 !!$ END DO 757 !!$ END DO 758 !!$ END DO 759 !!$ IF( lk_mpp ) CALL mppobc(tnbnd,jpind,jpinf,jpjnob+1,jpk*2*2,1,jpi, numout ) 760 !!$ IF( lk_mpp ) CALL mppobc(snbnd,jpind,jpinf,jpjnob+1,jpk*2*2,1,jpi, numout ) 761 !!$ 762 !!$ ! ... extremeties njn0,njn1 763 !!$ ii = jpind + 1 - nimpp 764 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 765 !!$ DO jk = 1, jpkm1 766 !!$ tnbnd(ii,jk,nibm,nitm) = tnbnd(ii+1,jk,nibm,nitm) 767 !!$ snbnd(ii,jk,nibm,nitm) = snbnd(ii+1,jk,nibm,nitm) 768 !!$ END DO 769 !!$ END IF 770 !!$ ii = jpinf + 1 - nimpp 771 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 772 !!$ DO jk = 1, jpkm1 773 !!$ tnbnd(ii,jk,nibm,nitm) = tnbnd(ii-1,jk,nibm,nitm) 774 !!$ snbnd(ii,jk,nibm,nitm) = snbnd(ii-1,jk,nibm,nitm) 775 !!$ END DO 776 !!$ END IF 777 !!$ 778 !!$ END IF ! End of array swap 779 !!$ 780 !!$ ! 2 - Calculation of radiation velocities 781 !!$ ! --------------------------------------- 782 !!$ 783 !!$ IF( kt >= nit000 +3 .OR. ln_rstart ) THEN 784 !!$ 785 !!$ ! 2.1 Calculate the normal velocity based on phase velocity u_cynbnd 786 !!$ ! ------------------------------------------------------------------- 787 !!$ ! 788 !!$ ! ji-row 789 !!$ ! | 790 !!$ ! nib -///u////// jpjnob + 1 791 !!$ ! /////|////// 792 !!$ ! nib -----f----- jpjnob 793 !!$ ! | 794 !!$ ! nibm-- u ---- jpjnob 795 !!$ ! | 796 !!$ ! nibm -----f----- jpjnob-1 797 !!$ ! | 798 !!$ ! nibm2-- u ---- jpjnob-1 799 !!$ ! | 800 !!$ ! nibm2 -----f----- jpjnob-2 801 !!$ ! | 802 !!$ ! ... radiative condition 803 !!$ ! ... jpjnob+1,(jpindp1, jpinfm1) 804 !!$ DO jj = fs_njn0+1, fs_njn1+1 ! Vector opt. 805 !!$ DO jk = 1, jpkm1 806 !!$ DO ji = 2, jpim1 807 !!$ ! ... 2* j-gradient of u (f-point i=nibm, time mean) 808 !!$ z2dx = ( unbnd(ji,jk,nibm ,nit) + unbnd(ji,jk,nibm ,nitm2) & 809 !!$ - 2.*unbnd(ji,jk,nibm2,nitm)) / e2f(ji,jj-2) 810 !!$ ! ... 2* i-gradient of u (u-point i=nibm, time nitm) 811 !!$ z2dy = ( unbnd(ji+1,jk,nibm,nitm) - unbnd(ji-1,jk,nibm,nitm) ) / e1u(ji,jj-1) 812 !!$ ! ... square of the norm of grad(v) 813 !!$ z4nor2 = z2dx * z2dx + z2dy * z2dy 814 !!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt 815 !!$ zdt = unbnd(ji,jk,nibm,nitm2) - unbnd(ji,jk,nibm,nit) 816 !!$ ! ... i-phase speed ratio (bounded by 1) and save the unbounded phase 817 !!$ ! velocity ratio no divided by e1f for the tracer radiation 818 !!$ IF( z4nor2 == 0.) THEN 819 !!$ z4nor2=.000001 820 !!$ END IF 821 !!$ z05cx = zdt * z2dx / z4nor2 822 !!$ u_cynbnd(ji,jk) = z05cx *unmsk(ji,jk) 823 !!$ END DO 824 !!$ END DO 825 !!$ END DO 826 !!$ IF( lk_mpp ) CALL mppobc(u_cynbnd,jpind,jpinf,jpjnob+1,jpk,1,jpi, numout ) 827 !!$ 828 !!$ ! ... extremeties njn0,njn1 829 !!$ ii = jpind + 1 - nimpp 830 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 831 !!$ DO jk = 1, jpkm1 832 !!$ u_cynbnd(ii,jk) = u_cynbnd(ii+1,jk) 833 !!$ END DO 834 !!$ END IF 835 !!$ ii = jpinf + 1 - nimpp 836 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 837 !!$ DO jk = 1, jpkm1 838 !!$ u_cynbnd(ii,jk) = u_cynbnd(ii-1,jk) 839 !!$ END DO 840 !!$ END IF 841 !!$ 842 !!$ ! 2.2 Calculate the normal velocity based on phase velocity v_cynbnd 843 !!$ ! ------------------------------------------------------------------ 844 !!$ ! 845 !!$ ! ji-row ji-row 846 !!$ ! | 847 !!$ ! /////|///////////////// 848 !!$ ! nib -----f----v----f---- jpjnob 849 !!$ ! | | 850 !!$ ! nib - u -- T -- u ---- jpjnob 851 !!$ ! | | 852 !!$ ! nibm -----f----v----f---- jpjnob-1 853 !!$ ! | | 854 !!$ ! nibm -- u -- T -- u --- jpjnob-1 855 !!$ ! | | 856 !!$ ! nibm2 -----f----v----f---- jpjnob-2 857 !!$ ! | | 858 !!$ ! ... Free surface formulation: 859 !!$ ! ... radiative conditions on the total part + relaxation toward climatology 860 !!$ ! ... jpjnob,(jpindp1, jpinfm1) 861 !!$ DO jj = fs_njn0, fs_njn1 ! Vector opt. 862 !!$ DO jk = 1, jpkm1 863 !!$ DO ji = 2, jpim1 864 !!$ ! ... 2* gradj(v) (T-point i=nibm, time mean) 865 !!$ ii = ji -1 + nimpp 866 !!$ z2dx = ( vnbnd(ji,jk,nibm ,nit) + vnbnd(ji,jk,nibm ,nitm2) & 867 !!$ - 2.*vnbnd(ji,jk,nibm2,nitm)) / e2t(ji,jj-1) 868 !!$ ! ... 2* gradi(v) (v-point i=nibm, time nitm) 869 !!$ z2dy = ( vnbnd(ji+1,jk,nibm,nitm) - vnbnd(ji-1,jk,nibm,nitm) ) / e1v(ji,jj-1) 870 !!$ ! ... square of the norm of grad(u) 871 !!$ z4nor2 = z2dx * z2dx + z2dy * z2dy 872 !!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt 873 !!$ zdt = vnbnd(ji,jk,nibm,nitm2) - vnbnd(ji,jk,nibm,nit) 874 !!$ ! ... j-phase speed ratio (bounded by 1) 875 !!$ IF( z4nor2 == 0. ) THEN 876 !!$ z4nor2=.00001 877 !!$ END IF 878 !!$ z05cx = zdt * z2dx / z4nor2 879 !!$ v_cynbnd(ji,jk)=z05cx *vnmsk(ji,jk) 880 !!$ END DO 881 !!$ END DO 882 !!$ END DO 883 !!$ 884 !!$ END IF 885 !!$ 886 !!$ END SUBROUTINE obc_rad_north 887 !!$ 888 !!$ 889 !!$ SUBROUTINE obc_rad_south ( kt ) 890 !!$ !!------------------------------------------------------------------------------ 891 !!$ !! *** SUBROUTINE obc_rad_south *** 892 !!$ !! 893 !!$ !! ** Purpose : 894 !!$ !! Perform swap of arrays to calculate radiative phase speeds at the open 895 !!$ !! south boundary and calculate those phase speeds if this OBC is not fixed. 896 !!$ !! In case of fixed OBC, this subrountine is not called. 897 !!$ !! 898 !!$ !! History : 899 !!$ !! ! 95-03 (J.-M. Molines) Original from SPEM 900 !!$ !! ! 97-07 (G. Madec, J.-M. Molines) additions 901 !!$ !! ! 97-12 (M. Imbard) Mpp adaptation 902 !!$ !! ! 00-06 (J.-M. Molines) 903 !!$ !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90 904 !!$ !!------------------------------------------------------------------------------ 905 !!$ !! * Arguments 906 !!$ INTEGER, INTENT( in ) :: kt 907 !!$ 908 !!$ !! * Local declarations 909 !!$ INTEGER :: ii 910 !!$ REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy 911 !!$ REAL(wp) :: zvcb, zvcbm, zvcbm2 912 !!$ !!------------------------------------------------------------------------------ 913 !!$ 914 !!$ ! 1. Swap arrays before calculating radiative velocities 915 !!$ ! ------------------------------------------------------ 916 !!$ 917 !!$ ! 1.1 zonal velocity 918 !!$ ! -------------------- 919 !!$ 920 !!$ IF( kt > nit000 .OR. ln_rstart ) THEN 921 !!$ 922 !!$ ! ... advance in time (time filter, array swap) 923 !!$ DO jk = 1, jpkm1 924 !!$ DO ji = 1, jpi 925 !!$ ! ... fields nitm2 <== nitm 926 !!$ usbnd(ji,jk,nib ,nitm2) = usbnd(ji,jk,nib ,nitm)*usmsk(ji,jk) 927 !!$ usbnd(ji,jk,nibm ,nitm2) = usbnd(ji,jk,nibm ,nitm)*usmsk(ji,jk) 928 !!$ usbnd(ji,jk,nibm2,nitm2) = usbnd(ji,jk,nibm2,nitm)*usmsk(ji,jk) 929 !!$ END DO 930 !!$ END DO 931 !!$ 932 !!$ DO jj = fs_njs0, fs_njs1 ! Vector opt. 933 !!$ DO jk = 1, jpkm1 934 !!$ DO ji = 1, jpi 935 !!$ usbnd(ji,jk,nib ,nitm) = usbnd(ji,jk,nib, nit)*usmsk(ji,jk) 936 !!$ usbnd(ji,jk,nibm ,nitm) = usbnd(ji,jk,nibm ,nit)*usmsk(ji,jk) 937 !!$ usbnd(ji,jk,nibm2,nitm) = usbnd(ji,jk,nibm2,nit)*usmsk(ji,jk) 938 !!$ ! ... fields nit <== now (kt+1) 939 !!$ usbnd(ji,jk,nib ,nit) = un(ji,jj ,jk)*usmsk(ji,jk) 940 !!$ usbnd(ji,jk,nibm ,nit) = un(ji,jj+1,jk)*usmsk(ji,jk) 941 !!$ usbnd(ji,jk,nibm2,nit) = un(ji,jj+2,jk)*usmsk(ji,jk) 942 !!$ END DO 943 !!$ END DO 944 !!$ END DO 945 !!$ IF( lk_mpp ) CALL mppobc(usbnd,jpisd,jpisf,jpjsob,jpk*3*3,1,jpi, numout ) 946 !!$ 947 !!$ ! ... extremeties njs0,njs1 948 !!$ ii = jpisd + 1 - nimpp 949 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 950 !!$ DO jk = 1, jpkm1 951 !!$ usbnd(ii,jk,nibm,nitm) = usbnd(ii+1,jk,nibm,nitm) 952 !!$ END DO 953 !!$ END IF 954 !!$ ii = jpisf + 1 - nimpp 955 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 956 !!$ DO jk = 1, jpkm1 957 !!$ usbnd(ii,jk,nibm,nitm) = usbnd(ii-1,jk,nibm,nitm) 958 !!$ END DO 959 !!$ END IF 960 !!$ 961 !!$ ! 1.2 normal velocity 962 !!$ ! ------------------- 963 !!$ 964 !!$ !.. advance in time (time filter, array swap) 965 !!$ DO jk = 1, jpkm1 966 !!$ DO ji = 1, jpi 967 !!$ ! ... fields nitm2 <== nitm 968 !!$ vsbnd(ji,jk,nib ,nitm2) = vsbnd(ji,jk,nib ,nitm)*vsmsk(ji,jk) 969 !!$ vsbnd(ji,jk,nibm ,nitm2) = vsbnd(ji,jk,nibm ,nitm)*vsmsk(ji,jk) 970 !!$ END DO 971 !!$ END DO 972 !!$ 973 !!$ DO jj = fs_njs0, fs_njs1 ! Vector opt. 974 !!$ DO jk = 1, jpkm1 975 !!$ DO ji = 1, jpi 976 !!$ vsbnd(ji,jk,nib ,nitm) = vsbnd(ji,jk,nib, nit)*vsmsk(ji,jk) 977 !!$ vsbnd(ji,jk,nibm ,nitm) = vsbnd(ji,jk,nibm ,nit)*vsmsk(ji,jk) 978 !!$ vsbnd(ji,jk,nibm2,nitm) = vsbnd(ji,jk,nibm2,nit)*vsmsk(ji,jk) 979 !!$ ! ... total or baroclinic velocity at b, bm and bm2 980 !!$ zvcb = vn (ji,jj,jk) 981 !!$ zvcbm = vn (ji,jj+1,jk) 982 !!$ zvcbm2 = vn (ji,jj+2,jk) 983 !!$ ! ... fields nit <== now (kt+1) 984 !!$ vsbnd(ji,jk,nib ,nit) = zvcb *vsmsk(ji,jk) 985 !!$ vsbnd(ji,jk,nibm ,nit) = zvcbm *vsmsk(ji,jk) 986 !!$ vsbnd(ji,jk,nibm2,nit) = zvcbm2 *vsmsk(ji,jk) 987 !!$ END DO 988 !!$ END DO 989 !!$ END DO 990 !!$ IF( lk_mpp ) CALL mppobc(vsbnd,jpisd,jpisf,jpjsob,jpk*3*3,1,jpi, numout ) 991 !!$ 992 !!$ ! ... extremeties njs0,njs1 993 !!$ ii = jpisd + 1 - nimpp 994 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 995 !!$ DO jk = 1, jpkm1 996 !!$ vsbnd(ii,jk,nibm,nitm) = vsbnd(ii+1,jk,nibm,nitm) 997 !!$ END DO 998 !!$ END IF 999 !!$ ii = jpisf + 1 - nimpp 1000 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 1001 !!$ DO jk = 1, jpkm1 1002 !!$ vsbnd(ii,jk,nibm,nitm) = vsbnd(ii-1,jk,nibm,nitm) 1003 !!$ END DO 1004 !!$ END IF 1005 !!$ 1006 !!$ ! 1.3 Temperature and salinity 1007 !!$ ! ---------------------------- 1008 !!$ 1009 !!$ ! ... advance in time (time filter, array swap) 1010 !!$ DO jk = 1, jpkm1 1011 !!$ DO ji = 1, jpi 1012 !!$ ! ... fields nitm <== nit plus time filter at the boundary 1013 !!$ tsbnd(ji,jk,nib,nitm) = tsbnd(ji,jk,nib,nit)*tsmsk(ji,jk) 1014 !!$ ssbnd(ji,jk,nib,nitm) = ssbnd(ji,jk,nib,nit)*tsmsk(ji,jk) 1015 !!$ END DO 1016 !!$ END DO 1017 !!$ 1018 !!$ DO jj = fs_njs0, fs_njs1 ! Vector opt. 1019 !!$ DO jk = 1, jpkm1 1020 !!$ DO ji = 1, jpi 1021 !!$ tsbnd(ji,jk,nibm ,nitm) = tsbnd(ji,jk,nibm ,nit)*tsmsk(ji,jk) 1022 !!$ ssbnd(ji,jk,nibm ,nitm) = ssbnd(ji,jk,nibm ,nit)*tsmsk(ji,jk) 1023 !!$ ! ... fields nit <== now (kt+1) 1024 !!$ tsbnd(ji,jk,nib ,nit) = tn(ji,jj ,jk)*tsmsk(ji,jk) 1025 !!$ tsbnd(ji,jk,nibm ,nit) = tn(ji,jj+1 ,jk)*tsmsk(ji,jk) 1026 !!$ ssbnd(ji,jk,nib ,nit) = sn(ji,jj ,jk)*tsmsk(ji,jk) 1027 !!$ ssbnd(ji,jk,nibm ,nit) = sn(ji,jj+1 ,jk)*tsmsk(ji,jk) 1028 !!$ END DO 1029 !!$ END DO 1030 !!$ END DO 1031 !!$ IF( lk_mpp ) CALL mppobc(tsbnd,jpisd,jpisf,jpjsob,jpk*2*2,1,jpi, numout ) 1032 !!$ IF( lk_mpp ) CALL mppobc(ssbnd,jpisd,jpisf,jpjsob,jpk*2*2,1,jpi, numout ) 1033 !!$ 1034 !!$ ! ... extremeties njs0,njs1 1035 !!$ ii = jpisd + 1 - nimpp 1036 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 1037 !!$ DO jk = 1, jpkm1 1038 !!$ tsbnd(ii,jk,nibm,nitm) = tsbnd(ii+1,jk,nibm,nitm) 1039 !!$ ssbnd(ii,jk,nibm,nitm) = ssbnd(ii+1,jk,nibm,nitm) 1040 !!$ END DO 1041 !!$ END IF 1042 !!$ ii = jpisf + 1 - nimpp 1043 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 1044 !!$ DO jk = 1, jpkm1 1045 !!$ tsbnd(ii,jk,nibm,nitm) = tsbnd(ii-1,jk,nibm,nitm) 1046 !!$ ssbnd(ii,jk,nibm,nitm) = ssbnd(ii-1,jk,nibm,nitm) 1047 !!$ END DO 1048 !!$ END IF 1049 !!$ 1050 !!$ END IF ! End of array swap 1051 !!$ 1052 !!$ ! 2 - Calculation of radiation velocities 1053 !!$ ! --------------------------------------- 1054 !!$ 1055 !!$ IF( kt >= nit000 +3 .OR. ln_rstart ) THEN 1056 !!$ 1057 !!$ ! 2.1 Calculate the normal velocity based on phase velocity u_cysbnd 1058 !!$ ! ------------------------------------------------------------------- 1059 !!$ ! 1060 !!$ ! ji-row 1061 !!$ ! | 1062 !!$ ! nibm2 -----f----- jpjsob +2 1063 !!$ ! | 1064 !!$ ! nibm2 -- u ----- jpjsob +2 1065 !!$ ! | 1066 !!$ ! nibm -----f----- jpjsob +1 1067 !!$ ! | 1068 !!$ ! nibm -- u ----- jpjsob +1 1069 !!$ ! | 1070 !!$ ! nib -----f----- jpjsob 1071 !!$ ! /////|////// 1072 !!$ ! nib ////u///// jpjsob 1073 !!$ ! 1074 !!$ ! ... radiative condition plus Raymond-Kuo 1075 !!$ ! ... jpjsob,(jpisdp1, jpisfm1) 1076 !!$ DO jj = fs_njs0, fs_njs1 ! Vector opt. 1077 !!$ DO jk = 1, jpkm1 1078 !!$ DO ji = 2, jpim1 1079 !!$ ! ... 2* j-gradient of u (f-point i=nibm, time mean) 1080 !!$ z2dx = (- usbnd(ji,jk,nibm ,nit) - usbnd(ji,jk,nibm ,nitm2) & 1081 !!$ + 2.*usbnd(ji,jk,nibm2,nitm) ) / e2f(ji,jj+1) 1082 !!$ ! ... 2* i-gradient of u (u-point i=nibm, time nitm) 1083 !!$ z2dy = ( usbnd(ji+1,jk,nibm,nitm) - usbnd(ji-1,jk,nibm,nitm) ) / e1u(ji, jj+1) 1084 !!$ ! ... square of the norm of grad(v) 1085 !!$ z4nor2 = z2dx * z2dx + z2dy * z2dy 1086 !!$ IF( z4nor2 == 0.) THEN 1087 !!$ z4nor2 = 0.000001 1088 !!$ END IF 1089 !!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt 1090 !!$ zdt = usbnd(ji,jk,nibm,nitm2) - usbnd(ji,jk,nibm,nit) 1091 !!$ ! ... i-phase speed ratio (bounded by -1) and save the unbounded phase 1092 !!$ ! velocity ratio no divided by e1f for the tracer radiation 1093 !!$ z05cx = zdt * z2dx / z4nor2 1094 !!$ u_cysbnd(ji,jk) = z05cx*usmsk(ji,jk) 1095 !!$ END DO 1096 !!$ END DO 1097 !!$ END DO 1098 !!$ IF( lk_mpp ) CALL mppobc(u_cysbnd,jpisd,jpisf,jpjsob,jpk,1,jpi, numout ) 1099 !!$ 1100 !!$ ! ... extremeties njs0,njs1 1101 !!$ ii = jpisd + 1 - nimpp 1102 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 1103 !!$ DO jk = 1, jpkm1 1104 !!$ u_cysbnd(ii,jk) = u_cysbnd(ii+1,jk) 1105 !!$ END DO 1106 !!$ END IF 1107 !!$ ii = jpisf + 1 - nimpp 1108 !!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN 1109 !!$ DO jk = 1, jpkm1 1110 !!$ u_cysbnd(ii,jk) = u_cysbnd(ii-1,jk) 1111 !!$ END DO 1112 !!$ END IF 1113 !!$ 1114 !!$ ! 2.2 Calculate the normal velocity based on phase velocity v_cysbnd 1115 !!$ ! ------------------------------------------------------------------- 1116 !!$ ! 1117 !!$ ! ji-row ji-row 1118 !!$ ! | | 1119 !!$ ! nibm2 -----f----v----f---- jpjsob+2 1120 !!$ ! | | 1121 !!$ ! nibm - u -- T -- u ---- jpjsob+2 1122 !!$ ! | | 1123 !!$ ! nibm -----f----v----f---- jpjsob+1 1124 !!$ ! | | 1125 !!$ ! nib -- u -- T -- u --- jpjsob+1 1126 !!$ ! | | 1127 !!$ ! nib -----f----v----f---- jpjsob 1128 !!$ ! ///////////////////// 1129 !!$ ! 1130 !!$ ! ... Free surface formulation: 1131 !!$ ! ... radiative conditions on the total part + relaxation toward climatology 1132 !!$ ! ... jpjsob,(jpisdp1,jpisfm1) 1133 !!$ DO jj = fs_njs0, fs_njs1 ! Vector opt. 1134 !!$ DO jk = 1, jpkm1 1135 !!$ DO ji = 2, jpim1 1136 !!$ ! ... 2* gradj(v) (T-point i=nibm, time mean) 1137 !!$ z2dx = ( - vsbnd(ji,jk,nibm ,nit) - vsbnd(ji,jk,nibm ,nitm2) & 1138 !!$ + 2.*vsbnd(ji,jk,nibm2,nitm) ) / e2t(ji,jj+1) 1139 !!$ ! ... 2* gradi(v) (v-point i=nibm, time nitm) 1140 !!$ z2dy = ( vsbnd(ji+1,jk,nibm,nitm) - vsbnd(ji-1,jk,nibm,nitm) ) / e1v(ji,jj+1) 1141 !!$ ! ... square of the norm of grad(u) 1142 !!$ z4nor2 = z2dx * z2dx + z2dy * z2dy 1143 !!$ IF( z4nor2 == 0.) THEN 1144 !!$ z4nor2 = 0.000001 1145 !!$ END IF 1146 !!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt 1147 !!$ zdt = vsbnd(ji,jk,nibm,nitm2) - vsbnd(ji,jk,nibm,nit) 1148 !!$ ! ... j-phase speed ratio (bounded by -1) 1149 !!$ z05cx = zdt * z2dx / z4nor2 1150 !!$ v_cysbnd(ji,jk)=z05cx*vsmsk(ji,jk) 1151 !!$ END DO 1152 !!$ END DO 1153 !!$ END DO 1154 !!$ 1155 !!$ ENDIF 1156 !!$ 1157 !!$ END SUBROUTINE obc_rad_south 1158 !!$ 1159 !!$#else 1160 1160 !!================================================================================= 1161 1161 !! *** MODULE obcrad ***
Note: See TracChangeset
for help on using the changeset viewer.
