- Timestamp:
- 2014-04-06T17:28:25+02:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r4596 r4616 60 60 # include "vectopt_loop_substitute.h90" 61 61 !!---------------------------------------------------------------------- 62 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)62 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 63 63 !! $Id$ 64 64 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 217 217 zcoef1 = zcoef0 * fse3w(ji,jj,1) 218 218 ! hydrostatic pressure gradient 219 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) /e1u(ji,jj)220 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) /e2v(ji,jj)219 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 220 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 221 221 ! add to the general momentum trend 222 222 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 234 234 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 235 235 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 236 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) /e1u(ji,jj)236 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 237 237 238 238 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 239 239 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 240 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) /e2v(ji,jj)240 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 241 241 ! add to the general momentum trend 242 242 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 284 284 zcoef1 = zcoef0 * fse3w(ji,jj,1) 285 285 ! hydrostatic pressure gradient 286 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) /e1u(ji,jj)287 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) /e2v(ji,jj)286 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 287 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 288 288 ! add to the general momentum trend 289 289 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 300 300 ! hydrostatic pressure gradient 301 301 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 302 & + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) &303 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) /e1u(ji,jj)302 & + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 303 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 304 304 305 305 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 306 & + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) &307 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) /e2v(ji,jj)306 & + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 307 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 308 308 ! add to the general momentum trend 309 309 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 315 315 316 316 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) 317 # if defined key_vectopt_loop318 jj = 1319 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)320 # else321 317 DO jj = 2, jpjm1 322 318 DO ji = 2, jpim1 323 # endif324 319 iku = mbku(ji,jj) 325 320 ikv = mbkv(ji,jj) … … 329 324 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value 330 325 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 331 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) /e1u(ji,jj)326 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 332 327 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 333 328 ENDIF … … 335 330 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value 336 331 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 337 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) /e2v(ji,jj)332 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 338 333 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 339 334 ENDIF 340 # if ! defined key_vectopt_loop 341 END DO 342 # endif 335 END DO 343 336 END DO 344 337 ! … … 392 385 DO ji = fs_2, fs_jpim1 ! vector opt. 393 386 ! hydrostatic pressure gradient along s-surfaces 394 zhpi(ji,jj,1) = zcoef0 /e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) &395 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) )396 zhpj(ji,jj,1) = zcoef0 /e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) &397 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) )387 zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 388 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 389 zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 390 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 398 391 ! s-coordinate pressure gradient correction 399 392 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 400 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) /e1u(ji,jj)393 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) * r1_e1u(ji,jj) 401 394 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 402 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) /e2v(ji,jj)395 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) * r1_e2v(ji,jj) 403 396 ! add to the general momentum trend 404 397 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 406 399 END DO 407 400 END DO 401 402 !!gm 403 !!gm Idea of optimization here : only divide by e1u and e2v when apply to ua, va .... 404 !!gm 408 405 409 406 ! interior value (2=<jk=<jpkm1) … … 412 409 DO ji = fs_2, fs_jpim1 ! vector opt. 413 410 ! hydrostatic pressure gradient along s-surfaces 414 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 /e1u(ji,jj) &411 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 415 412 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 416 413 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 417 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 /e2v(ji,jj) &414 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 418 415 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 419 416 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 420 417 ! s-coordinate pressure gradient correction 421 418 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 422 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) /e1u(ji,jj)419 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) * r1_e1u(ji,jj) 423 420 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 424 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) /e2v(ji,jj)421 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) * r1_e2v(ji,jj) 425 422 ! add to the general momentum trend 426 423 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 625 622 DO jj = 2, jpjm1 626 623 DO ji = fs_2, fs_jpim1 ! vector opt. 627 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) /e1u(ji,jj)628 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) /e2v(ji,jj)624 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 625 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 629 626 ! add to the general momentum trend 630 627 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 642 639 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 643 640 & + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) & 644 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) /e1u(ji,jj)641 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj) 645 642 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 646 643 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 647 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) /e2v(ji,jj)644 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 648 645 ! add to the general momentum trend 649 646 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 853 850 ! update the momentum trends in u direction 854 851 855 zdpdx1 = zcoef0 /e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))852 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 856 853 IF( lk_vvl ) THEN 857 zdpdx2 = zcoef0 /e1u(ji,jj) * &854 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 858 855 ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 859 856 ELSE 860 zdpdx2 = zcoef0 /e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)857 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 861 858 ENDIF 862 859 863 860 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 864 861 & umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 862 !!gm above it is stupid: umask is equal to the product of the 2 tmask that follows... 865 863 ENDIF 866 864 … … 910 908 ! update the momentum trends in v direction 911 909 912 zdpdy1 = zcoef0 /e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))910 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 913 911 IF( lk_vvl ) THEN 914 zdpdy2 = zcoef0 /e2v(ji,jj) * &912 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 915 913 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 916 914 ELSE … … 920 918 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 921 919 & vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 920 !!gm above it is stupid: vmask is equal to the product of the 2 tmask that follows... 922 921 ENDIF 923 922 … … 932 931 END SUBROUTINE hpg_prj 933 932 933 934 934 SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 935 935 !!---------------------------------------------------------------------- … … 940 940 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 941 941 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 942 !! 943 !!---------------------------------------------------------------------- 944 IMPLICIT NONE 942 !!---------------------------------------------------------------------- 945 943 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate 946 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 947 ! the interpoated function 944 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of the interpolated function 948 945 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 949 ! 2: Linear 950 951 ! Local Variables 946 ! ! 2: Linear 947 ! 952 948 INTEGER :: ji, jj, jk ! dummy loop indices 953 949 INTEGER :: jpi, jpj, jpkm1 … … 1038 1034 CALL ctl_stop( 'invalid polynomial type in cspline' ) 1039 1035 ENDIF 1040 1041 1036 ! 1042 1037 END SUBROUTINE cspline 1043 1038 … … 1069 1064 END FUNCTION interp1 1070 1065 1066 1071 1067 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1072 1068 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
