Changeset 3598 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
- Timestamp:
- 2012-11-19T14:35:09+01:00 (11 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 (modified) (56 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r3434 r3598 11 11 !! 8.5 ! 2002-08 (A. Bozec) hpg_zps: Original code 12 12 !! NEMO 1.0 ! 2005-10 (A. Beckmann, B.W. An) various s-coordinate options 13 !! ! Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot 13 !! ! Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot 14 14 !! - ! 2005-11 (G. Madec) style & small optimisation 15 15 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase … … 31 31 USE dom_oce ! ocean space and time domain 32 32 USE phycst ! physical constants 33 USE trdmod ! ocean dynamics trends 33 USE trdmod ! ocean dynamics trends 34 34 USE trdmod_oce ! ocean variables trends 35 35 USE in_out_manager ! I/O manager 36 36 USE prtctl ! Print control 37 USE lbclnk ! lateral boundary condition 37 USE lbclnk ! lateral boundary condition 38 38 USE lib_mpp ! MPP library 39 39 USE wrk_nemo ! Memory Allocation … … 46 46 PUBLIC dyn_hpg_init ! routine called by opa module 47 47 48 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 48 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 49 49 LOGICAL , PUBLIC :: ln_hpg_zco = .TRUE. !: z-coordinate - full steps 50 50 LOGICAL , PUBLIC :: ln_hpg_zps = .FALSE. !: z-coordinate - partial steps (interpolation) … … 54 54 LOGICAL , PUBLIC :: ln_dynhpg_imp = .FALSE. !: semi-implicite hpg flag 55 55 56 INTEGER :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags)56 INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 57 57 58 58 !! * Substitutions … … 70 70 !! *** ROUTINE dyn_hpg *** 71 71 !! 72 !! ** Method : Call the hydrostatic pressure gradient routine 72 !! ** Method : Call the hydrostatic pressure gradient routine 73 73 !! using the scheme defined in the namelist 74 !! 74 !! 75 75 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 76 76 !! - Save the trend (l_trddyn=T) … … 84 84 IF( l_trddyn ) THEN ! Temporary saving of ua and va trends (l_trddyn) 85 85 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 86 ztrdu(:,:,:) = ua(:,:,:) 87 ztrdv(:,:,:) = va(:,:,:) 88 ENDIF 86 ztrdu(:,:,:) = ua(:,:,:) 87 ztrdv(:,:,:) = va(:,:,:) 88 ENDIF 89 89 ! 90 90 SELECT CASE ( nhpg ) ! Hydrostatic pressure gradient computation … … 101 101 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt ) 102 102 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 103 ENDIF 103 ENDIF 104 104 ! 105 105 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg - Ua: ', mask1=umask, & … … 161 161 ! 162 162 ! ! Consistency check 163 ioptio = 0 163 ioptio = 0 164 164 IF( ln_hpg_zco ) ioptio = ioptio + 1 165 165 IF( ln_hpg_zps ) ioptio = ioptio + 1 … … 185 185 !! ua = ua - 1/e1u * zhpi 186 186 !! va = va - 1/e2v * zhpj 187 !! 187 !! 188 188 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 189 189 !!---------------------------------------------------------------------- … … 192 192 INTEGER :: ji, jj, jk ! dummy loop indices 193 193 REAL(wp) :: zcoef0, zcoef1 ! temporary scalars 194 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 195 !!---------------------------------------------------------------------- 196 ! 194 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 195 !!---------------------------------------------------------------------- 196 ! 197 197 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 198 198 ! … … 202 202 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case ' 203 203 ENDIF 204 205 zcoef0 = - grav * 0.5_wp ! Local constant initialization 204 205 zcoef0 = - grav * 0.5_wp ! Local constant initialization 206 206 207 207 ! Surface value … … 247 247 !!--------------------------------------------------------------------- 248 248 !! *** ROUTINE hpg_zps *** 249 !! 249 !! 250 250 !! ** Method : z-coordinate plus partial steps case. blahblah... 251 !! 251 !! 252 252 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 253 !!---------------------------------------------------------------------- 253 !!---------------------------------------------------------------------- 254 254 INTEGER, INTENT(in) :: kt ! ocean time-step index 255 255 !! … … 257 257 INTEGER :: iku, ikv ! temporary integers 258 258 REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars 259 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 259 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 260 260 !!---------------------------------------------------------------------- 261 261 ! … … 363 363 INTEGER :: ji, jj, jk ! dummy loop indices 364 364 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 365 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 365 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 366 366 !!---------------------------------------------------------------------- 367 367 ! … … 383 383 ! Surface value 384 384 DO jj = 2, jpjm1 385 DO ji = fs_2, fs_jpim1 ! vector opt. 385 DO ji = fs_2, fs_jpim1 ! vector opt. 386 386 ! hydrostatic pressure gradient along s-surfaces 387 387 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & … … 397 397 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 398 398 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 399 END DO 400 END DO 401 399 END DO 400 END DO 401 402 402 ! interior value (2=<jk=<jpkm1) 403 DO jk = 2, jpkm1 404 DO jj = 2, jpjm1 405 DO ji = fs_2, fs_jpim1 ! vector opt. 403 DO jk = 2, jpkm1 404 DO jj = 2, jpjm1 405 DO ji = fs_2, fs_jpim1 ! vector opt. 406 406 ! hydrostatic pressure gradient along s-surfaces 407 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 408 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 407 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 408 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 409 409 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 410 410 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & … … 432 432 !! 433 433 !! ** Method : Density Jacobian with Cubic polynomial scheme 434 !! 434 !! 435 435 !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003 436 436 !!---------------------------------------------------------------------- … … 441 441 REAL(wp) :: z1_10, cffu, cffx ! " " 442 442 REAL(wp) :: z1_12, cffv, cffy ! " " 443 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 443 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 444 444 REAL(wp), POINTER, DIMENSION(:,:,:) :: dzx, dzy, dzz, dzu, dzv, dzw 445 445 REAL(wp), POINTER, DIMENSION(:,:,:) :: drhox, drhoy, drhoz, drhou, drhov, drhow … … 447 447 !!---------------------------------------------------------------------- 448 448 ! 449 CALL wrk_alloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw ) 450 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 451 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 449 CALL wrk_alloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw ) 450 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 451 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 452 452 ! 453 453 … … 497 497 cffu = 2._wp * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk ) 498 498 cffx = 2._wp * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk ) 499 499 500 500 cffv = 2._wp * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk ) 501 501 cffy = 2._wp * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk ) … … 568 568 & + 0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) ) & 569 569 & * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) ) & 570 & / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) ) 570 & / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) ) 571 571 END DO 572 572 END DO … … 631 631 ! ---------------- 632 632 DO jk = 2, jpkm1 633 DO jj = 2, jpjm1 633 DO jj = 2, jpjm1 634 634 DO ji = fs_2, fs_jpim1 ! vector opt. 635 635 ! hydrostatic pressure gradient along s-surfaces … … 647 647 END DO 648 648 ! 649 CALL wrk_dealloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw ) 650 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 651 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 649 CALL wrk_dealloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw ) 650 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 651 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 652 652 ! 653 653 END SUBROUTINE hpg_djc … … 676 676 INTEGER :: jk1, jis, jid, jjs, jjd 677 677 REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 678 REAL(wp) :: zrhdt1 678 REAL(wp) :: zrhdt1 679 679 REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 680 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdept, zrhh 680 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdept, zrhh 681 681 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 682 682 !!---------------------------------------------------------------------- 683 683 ! 684 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 685 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 684 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 685 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 686 686 ! 687 687 IF( kt == nit000 ) THEN … … 693 693 !!---------------------------------------------------------------------- 694 694 ! Local constant initialization 695 zcoef0 = - grav 695 zcoef0 = - grav 696 696 znad = 0.0_wp 697 697 IF( lk_vvl ) znad = 1._wp … … 700 700 zhpi(:,:,:) = 0._wp 701 701 zrhh(:,:,:) = rhd(:,:,:) 702 702 703 703 ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate 704 704 DO jj = 1, jpj 705 DO ji = 1, jpi 705 DO ji = 1, jpi 706 706 jk = mbathy(ji,jj) 707 707 IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp … … 711 711 zrhh(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1), & 712 712 fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 713 END DO 713 END DO 714 714 ENDIF 715 715 END DO … … 728 728 xsp(:,:,:) = zdept(:,:,:) 729 729 730 ! Construct the vertical density profile with the 730 ! Construct the vertical density profile with the 731 731 ! constrained cubic spline interpolation 732 732 ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 733 CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type) 733 CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type) 734 734 735 735 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 736 736 DO jj = 2, jpj 737 DO ji = 2, jpi 737 DO ji = 2, jpi 738 738 zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), & 739 739 bsp(ji,jj,1), csp(ji,jj,1), & … … 741 741 742 742 ! assuming linear profile across the top half surface layer 743 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * zrhdt1 743 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * zrhdt1 744 744 END DO 745 745 END DO 746 746 747 747 ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)" 748 DO jk = 2, jpkm1 749 DO jj = 2, jpj 748 DO jk = 2, jpkm1 749 DO jj = 2, jpj 750 750 DO ji = 2, jpi 751 751 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & … … 758 758 759 759 ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 760 DO jj = 2, jpjm1 761 DO ji = 2, jpim1 760 DO jj = 2, jpjm1 761 DO ji = 2, jpim1 762 762 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad) 763 763 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad) … … 765 765 END DO 766 766 767 DO jk = 2, jpkm1 768 DO jj = 2, jpjm1 769 DO ji = 2, jpim1 767 DO jk = 2, jpkm1 768 DO jj = 2, jpjm1 769 DO ji = 2, jpim1 770 770 zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk) 771 771 zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk) … … 773 773 END DO 774 774 END DO 775 776 DO jk = 1, jpkm1 777 DO jj = 2, jpjm1 778 DO ji = 2, jpim1 775 776 DO jk = 1, jpkm1 777 DO jj = 2, jpjm1 778 DO ji = 2, jpim1 779 779 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk) 780 780 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk) … … 795 795 796 796 797 DO jk = 1, jpkm1 798 DO jj = 2, jpjm1 799 DO ji = 2, jpim1 797 DO jk = 1, jpkm1 798 DO jj = 2, jpjm1 799 DO ji = 2, jpim1 800 800 zpwes = 0._wp; zpwed = 0._wp 801 801 zpnss = 0._wp; zpnsd = 0._wp … … 812 812 813 813 ! integrate the pressure on the shallow side 814 jk1 = jk 814 jk1 = jk 815 815 DO WHILE ( -zdept(jis,jj,jk1) > zuijk ) 816 816 IF( jk1 == mbku(ji,jj) ) THEN … … 819 819 ENDIF 820 820 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk) 821 zpwes = zpwes + & 821 zpwes = zpwes + & 822 822 integ_spline(zdept(jis,jj,jk1), zdeps, & 823 823 asp(jis,jj,jk1), bsp(jis,jj,jk1), & … … 825 825 jk1 = jk1 + 1 826 826 END DO 827 827 828 828 ! integrate the pressure on the deep side 829 jk1 = jk 829 jk1 = jk 830 830 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 831 831 IF( jk1 == 1 ) THEN … … 838 838 ENDIF 839 839 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk) 840 zpwed = zpwed + & 840 zpwed = zpwed + & 841 841 integ_spline(zdeps, zdept(jid,jj,jk1), & 842 842 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & … … 844 844 jk1 = jk1 - 1 845 845 END DO 846 846 847 847 ! update the momentum trends in u direction 848 848 849 849 zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 850 850 IF( lk_vvl ) THEN 851 zdpdx2 = zcoef0 / e1u(ji,jj) * & 852 ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 851 zdpdx2 = zcoef0 / e1u(ji,jj) * & 852 ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 853 853 ELSE 854 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 854 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 855 855 ENDIF 856 856 … … 858 858 & umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 859 859 ENDIF 860 860 861 861 !!!!! for v equation 862 862 IF( jk <= mbkv(ji,jj) ) THEN … … 868 868 869 869 ! integrate the pressure on the shallow side 870 jk1 = jk 870 jk1 = jk 871 871 DO WHILE ( -zdept(ji,jjs,jk1) > zvijk ) 872 872 IF( jk1 == mbkv(ji,jj) ) THEN … … 875 875 ENDIF 876 876 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk) 877 zpnss = zpnss + & 877 zpnss = zpnss + & 878 878 integ_spline(zdept(ji,jjs,jk1), zdeps, & 879 879 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & … … 881 881 jk1 = jk1 + 1 882 882 END DO 883 883 884 884 ! integrate the pressure on the deep side 885 jk1 = jk 885 jk1 = jk 886 886 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 887 887 IF( jk1 == 1 ) THEN … … 894 894 ENDIF 895 895 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk) 896 zpnsd = zpnsd + & 896 zpnsd = zpnsd + & 897 897 integ_spline(zdeps, zdept(ji,jjd,jk1), & 898 898 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & … … 900 900 jk1 = jk1 - 1 901 901 END DO 902 902 903 903 904 904 ! update the momentum trends in v direction … … 907 907 IF( lk_vvl ) THEN 908 908 zdpdy2 = zcoef0 / e2v(ji,jj) * & 909 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 909 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 910 910 ELSE 911 zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 911 zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 912 912 ENDIF 913 913 … … 916 916 ENDIF 917 917 918 918 919 919 END DO 920 920 END DO 921 921 END DO 922 922 ! 923 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 924 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 923 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 924 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 925 925 ! 926 926 END SUBROUTINE hpg_prj … … 929 929 !!---------------------------------------------------------------------- 930 930 !! *** ROUTINE cspline *** 931 !! 931 !! 932 932 !! ** Purpose : constrained cubic spline interpolation 933 !! 934 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 933 !! 934 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 935 935 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 936 936 !! … … 938 938 IMPLICIT NONE 939 939 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate 940 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 940 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 941 941 ! the interpoated function 942 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 942 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 943 943 ! 2: Linear 944 944 945 ! Local Variables 945 ! Local Variables 946 946 INTEGER :: ji, jj, jk ! dummy loop indices 947 947 INTEGER :: jpi, jpj, jpkm1 … … 955 955 jpkm1 = size(fsp,3) - 1 956 956 957 957 958 958 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 959 959 DO ji = 1, jpi 960 960 DO jj = 1, jpj 961 !!Fritsch&Butland's method, 1984 (preferred, but more computation) 961 !!Fritsch&Butland's method, 1984 (preferred, but more computation) 962 962 ! DO jk = 2, jpkm1-1 963 ! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 964 ! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 963 ! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 964 ! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 965 965 ! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 966 966 ! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 967 967 ! 968 968 ! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 969 ! 969 ! 970 970 ! IF(zdf1 * zdf2 <= 0._wp) THEN 971 971 ! zdf(jk) = 0._wp … … 974 974 ! ENDIF 975 975 ! END DO 976 976 977 977 !!Simply geometric average 978 978 DO jk = 2, jpkm1-1 979 979 zdf1 = (fsp(ji,jj,jk) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk) - xsp(ji,jj,jk-1)) 980 980 zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk)) 981 981 982 982 IF(zdf1 * zdf2 <= 0._wp) THEN 983 983 zdf(jk) = 0._wp … … 986 986 ENDIF 987 987 END DO 988 988 989 989 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 990 990 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) … … 992 992 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 993 993 & 0.5_wp * zdf(jpkm1 - 1) 994 994 995 995 DO jk = 1, jpkm1 - 1 996 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 996 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 997 997 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 998 998 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 999 zddf1 = -2._wp * ztmp1 + ztmp2 999 zddf1 = -2._wp * ztmp1 + ztmp2 1000 1000 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1001 zddf2 = 2._wp * ztmp1 - ztmp2 1002 1001 zddf2 = 2._wp * ztmp1 - ztmp2 1002 1003 1003 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1004 1004 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1005 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1005 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1006 1006 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1007 1007 & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - & … … 1013 1013 END DO 1014 1014 END DO 1015 1015 1016 1016 ELSE IF (polynomial_type == 2) THEN ! Linear 1017 1017 DO ji = 1, jpi 1018 1018 DO jj = 1, jpj 1019 1019 DO jk = 1, jpkm1-1 1020 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1020 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1021 1021 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1022 1022 1023 1023 dsp(ji,jj,jk) = 0._wp 1024 1024 csp(ji,jj,jk) = 0._wp … … 1033 1033 ENDIF 1034 1034 1035 1035 1036 1036 END SUBROUTINE cspline 1037 1037 1038 1038 1039 FUNCTION interp1(x, xl, xr, fl, fr) RESULT(f) 1039 FUNCTION interp1(x, xl, xr, fl, fr) RESULT(f) 1040 1040 !!---------------------------------------------------------------------- 1041 1041 !! *** ROUTINE interp1 *** 1042 !! 1042 !! 1043 1043 !! ** Purpose : 1-d linear interpolation 1044 !! 1045 !! ** Method : 1044 !! 1045 !! ** Method : 1046 1046 !! interpolation is straight forward 1047 !! extrapolation is also permitted (no value limit) 1047 !! extrapolation is also permitted (no value limit) 1048 1048 !! 1049 1049 !!---------------------------------------------------------------------- 1050 1050 IMPLICIT NONE 1051 REAL(wp), INTENT(in) :: x, xl, xr, fl, fr 1051 REAL(wp), INTENT(in) :: x, xl, xr, fl, fr 1052 1052 REAL(wp) :: f ! result of the interpolation (extrapolation) 1053 1053 REAL(wp) :: zdeltx … … 1060 1060 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 1061 1061 ENDIF 1062 1062 1063 1063 END FUNCTION interp1 1064 1064 1065 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1065 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1066 1066 !!---------------------------------------------------------------------- 1067 1067 !! *** ROUTINE interp1 *** 1068 !! 1068 !! 1069 1069 !! ** Purpose : 1-d constrained cubic spline interpolation 1070 !! 1070 !! 1071 1071 !! ** Method : cubic spline interpolation 1072 1072 !! 1073 1073 !!---------------------------------------------------------------------- 1074 1074 IMPLICIT NONE 1075 REAL(wp), INTENT(in) :: x, a, b, c, d 1075 REAL(wp), INTENT(in) :: x, a, b, c, d 1076 1076 REAL(wp) :: f ! value from the interpolation 1077 1077 !!---------------------------------------------------------------------- 1078 1078 1079 f = a + x* ( b + x * ( c + d * x ) ) 1079 f = a + x* ( b + x * ( c + d * x ) ) 1080 1080 1081 1081 END FUNCTION interp2 1082 1082 1083 1083 1084 FUNCTION interp3(x, a, b, c, d) RESULT(f) 1084 FUNCTION interp3(x, a, b, c, d) RESULT(f) 1085 1085 !!---------------------------------------------------------------------- 1086 1086 !! *** ROUTINE interp1 *** 1087 !! 1087 !! 1088 1088 !! ** Purpose : Calculate the first order of deriavtive of 1089 1089 !! a cubic spline function y=a+b*x+c*x^2+d*x^3 1090 !! 1090 !! 1091 1091 !! ** Method : f=dy/dx=b+2*c*x+3*d*x^2 1092 1092 !! 1093 1093 !!---------------------------------------------------------------------- 1094 1094 IMPLICIT NONE 1095 REAL(wp), INTENT(in) :: x, a, b, c, d 1095 REAL(wp), INTENT(in) :: x, a, b, c, d 1096 1096 REAL(wp) :: f ! value from the interpolation 1097 1097 !!---------------------------------------------------------------------- … … 1101 1101 END FUNCTION interp3 1102 1102 1103 1104 FUNCTION integ_spline(xl, xr, a, b, c, d) RESULT(f) 1103 1104 FUNCTION integ_spline(xl, xr, a, b, c, d) RESULT(f) 1105 1105 !!---------------------------------------------------------------------- 1106 1106 !! *** ROUTINE interp1 *** 1107 !! 1107 !! 1108 1108 !! ** Purpose : 1-d constrained cubic spline integration 1109 !! 1110 !! ** Method : integrate polynomial a+bx+cx^2+dx^3 from xl to xr 1109 !! 1110 !! ** Method : integrate polynomial a+bx+cx^2+dx^3 from xl to xr 1111 1111 !! 1112 1112 !!---------------------------------------------------------------------- 1113 1113 IMPLICIT NONE 1114 REAL(wp), INTENT(in) :: xl, xr, a, b, c, d 1115 REAL(wp) :: za1, za2, za3 1114 REAL(wp), INTENT(in) :: xl, xr, a, b, c, d 1115 REAL(wp) :: za1, za2, za3 1116 1116 REAL(wp) :: f ! integration result 1117 1117 !!---------------------------------------------------------------------- 1118 1118 1119 za1 = 0.5_wp * b 1120 za2 = c / 3.0_wp 1121 za3 = 0.25_wp * d 1119 za1 = 0.5_wp * b 1120 za2 = c / 3.0_wp 1121 za3 = 0.25_wp * d 1122 1122 1123 1123 f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - &
Note: See TracChangeset
for help on using the changeset viewer.
