- Timestamp:
- 2012-11-19T14:35:09+01:00 (10 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 5 edited
-
DYN/dynhpg.F90 (modified) (56 diffs)
-
LBC/lib_mpp.F90 (modified) (134 diffs)
-
SBC/sbcrnf.F90 (modified) (24 diffs)
-
TRA/trabbl.F90 (modified) (40 diffs)
-
ZDF/zdfbfr.F90 (modified) (8 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 ) ) ) - & -
trunk/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
r3435 r3598 17 17 !! - ! 2008 (R. Benshila) add mpp_ini_ice 18 18 !! 3.2 ! 2009 (R. Benshila) SHMEM suppression, north fold in lbc_nfd 19 !! 3.2 ! 2009 (O. Marti) add mpp_ini_znl 19 !! 3.2 ! 2009 (O. Marti) add mpp_ini_znl 20 20 !! 4.0 ! 2011 (G. Madec) move ctl_ routines from in_out_manager 21 21 !!---------------------------------------------------------------------- … … 27 27 !! get_unit : give the index of an unused logical unit 28 28 !!---------------------------------------------------------------------- 29 #if defined key_mpp_mpi 29 #if defined key_mpp_mpi 30 30 !!---------------------------------------------------------------------- 31 31 !! 'key_mpp_mpi' MPI massively parallel processing library … … 52 52 !! mpp_lbc_north_e : variant of mpp_lbc_north for extra outer halo 53 53 !!---------------------------------------------------------------------- 54 USE dom_oce ! ocean space and time domain 54 USE dom_oce ! ocean space and time domain 55 55 USE lbcnfd ! north fold treatment 56 56 USE in_out_manager ! I/O manager … … 58 58 IMPLICIT NONE 59 59 PRIVATE 60 60 61 61 PUBLIC ctl_stop, ctl_warn, get_unit, ctl_opn 62 62 PUBLIC mynode, mppstop, mppsync, mpp_comm_free … … 67 67 PUBLIC mppobc, mpp_ini_ice, mpp_ini_znl 68 68 PUBLIC mppsize 69 PUBLIC lib_mpp_alloc ! Called in nemogcm.F90 69 PUBLIC lib_mpp_alloc ! Called in nemogcm.F90 70 PUBLIC mppsend, mpprecv ! (PUBLIC for TAM) 70 71 71 72 !! * Interfaces … … 84 85 END INTERFACE 85 86 INTERFACE mpp_lbc_north 86 MODULE PROCEDURE mpp_lbc_north_3d, mpp_lbc_north_2d 87 MODULE PROCEDURE mpp_lbc_north_3d, mpp_lbc_north_2d 87 88 END INTERFACE 88 89 INTERFACE mpp_minloc … … 92 93 MODULE PROCEDURE mpp_maxloc2d ,mpp_maxloc3d 93 94 END INTERFACE 94 95 95 96 !! ========================= !! 96 97 !! MPI variable definition !! … … 99 100 INCLUDE 'mpif.h' 100 101 !$AGRIF_END_DO_NOT_TREAT 101 102 102 103 LOGICAL, PUBLIC, PARAMETER :: lk_mpp = .TRUE. !: mpp flag 103 104 104 105 INTEGER, PARAMETER :: nprocmax = 2**10 ! maximun dimension (required to be a power of 2) 105 106 106 107 INTEGER :: mppsize ! number of process 107 108 INTEGER :: mpprank ! process number [ 0 - size-1 ] … … 126 127 INTEGER :: ndim_rank_znl ! number of processors on the same zonal average 127 128 INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: nrank_znl ! dimension ndim_rank_znl, number of the procs into the same znl domain 128 129 129 130 ! North fold condition in mpp_mpi with jpni > 1 130 131 INTEGER :: ngrp_world ! group ID for the world processors … … 140 141 CHARACTER(len=1) :: cn_mpi_send = 'S' ! type od mpi send/recieve (S=standard, B=bsend, I=isend) 141 142 LOGICAL :: l_isend = .FALSE. ! isend use indicator (T if cn_mpi_send='I') 142 INTEGER :: nn_buffer = 0 ! size of the buffer in case of mpi_bsend 143 143 INTEGER :: nn_buffer = 0 ! size of the buffer in case of mpi_bsend 144 144 145 REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: tampon ! buffer in case of bsend 145 146 … … 173 174 ! North fold arrays used to minimise the use of allgather operations. Set in nemo_northcomms (nemogcm) so need to be public 174 175 INTEGER, PUBLIC, PARAMETER :: jpmaxngh = 8 ! Assumed maximum number of active neighbours 175 INTEGER, PUBLIC, PARAMETER :: jptyps = 5 ! Number of different neighbour lists to be used for northfold exchanges 176 INTEGER, PUBLIC, PARAMETER :: jptyps = 5 ! Number of different neighbour lists to be used for northfold exchanges 176 177 INTEGER, PUBLIC, DIMENSION (jpmaxngh,jptyps) :: isendto 177 178 INTEGER, PUBLIC, DIMENSION (jptyps) :: nsndto … … 229 230 !!---------------------------------------------------------------------- 230 231 !! *** routine mynode *** 231 !! 232 !! 232 233 !! ** Purpose : Find processor unit 233 234 !!---------------------------------------------------------------------- 234 CHARACTER(len=*),DIMENSION(:), INTENT( out) :: ldtxt 235 INTEGER , INTENT(in ) :: kumnam ! namelist logical unit 236 INTEGER , INTENT(inout) :: kstop ! stop indicator 235 CHARACTER(len=*),DIMENSION(:), INTENT( out) :: ldtxt 236 INTEGER , INTENT(in ) :: kumnam ! namelist logical unit 237 INTEGER , INTENT(inout) :: kstop ! stop indicator 237 238 INTEGER, OPTIONAL , INTENT(in ) :: localComm 238 239 ! … … 258 259 #if defined key_agrif 259 260 IF( .NOT. Agrif_Root() ) THEN 260 jpni = Agrif_Parent(jpni ) 261 jpni = Agrif_Parent(jpni ) 261 262 jpnj = Agrif_Parent(jpnj ) 262 263 jpnij = Agrif_Parent(jpnij) … … 282 283 CALL mpi_initialized ( mpi_was_called, code ) 283 284 IF( code /= MPI_SUCCESS ) THEN 284 DO ji = 1, SIZE(ldtxt) 285 DO ji = 1, SIZE(ldtxt) 285 286 IF( TRIM(ldtxt(ji)) /= '' ) WRITE(*,*) ldtxt(ji) ! control print of mynode 286 END DO 287 END DO 287 288 WRITE(*, cform_err) 288 289 WRITE(*, *) 'lib_mpp: Error in routine mpi_initialized' … … 297 298 CASE ( 'B' ) ! Buffer mpi send (blocking) 298 299 WRITE(ldtxt(ii),*) ' Buffer blocking mpi send (bsend)' ; ii = ii + 1 299 IF( Agrif_Root() ) CALL mpi_init_opa( ldtxt, ii, ierr ) 300 IF( Agrif_Root() ) CALL mpi_init_opa( ldtxt, ii, ierr ) 300 301 CASE ( 'I' ) ! Immediate mpi send (non-blocking send) 301 302 WRITE(ldtxt(ii),*) ' Immediate non-blocking send (isend)' ; ii = ii + 1 … … 330 331 ENDIF 331 332 332 IF( PRESENT(localComm) ) THEN 333 IF( PRESENT(localComm) ) THEN 333 334 IF( Agrif_Root() ) THEN 334 335 mpi_comm_opa = localComm … … 337 338 CALL mpi_comm_dup( mpi_comm_world, mpi_comm_opa, code) 338 339 IF( code /= MPI_SUCCESS ) THEN 339 DO ji = 1, SIZE(ldtxt) 340 DO ji = 1, SIZE(ldtxt) 340 341 IF( TRIM(ldtxt(ji)) /= '' ) WRITE(*,*) ldtxt(ji) ! control print of mynode 341 342 END DO … … 344 345 CALL mpi_abort( mpi_comm_world, code, ierr ) 345 346 ENDIF 346 ENDIF 347 ENDIF 347 348 348 349 CALL mpi_comm_rank( mpi_comm_opa, mpprank, ierr ) 349 350 CALL mpi_comm_size( mpi_comm_opa, mppsize, ierr ) 350 351 mynode = mpprank 351 ! 352 ! 352 353 CALL MPI_OP_CREATE(DDPDD_MPI, .TRUE., MPI_SUMDD, ierr) 353 354 ! … … 361 362 !! ** Purpose : Message passing manadgement 362 363 !! 363 !! ** Method : Use mppsend and mpprecv function for passing mask 364 !! ** Method : Use mppsend and mpprecv function for passing mask 364 365 !! between processors following neighboring subdomains. 365 366 !! domain parameters … … 368 369 !! nbondi : mark for "east-west local boundary" 369 370 !! nbondj : mark for "north-south local boundary" 370 !! noea : number for local neighboring processors 371 !! noea : number for local neighboring processors 371 372 !! nowe : number for local neighboring processors 372 373 !! noso : number for local neighboring processors … … 381 382 REAL(wp) , INTENT(in ) :: psgn ! =-1 the sign change across the north fold boundary 382 383 ! ! = 1. , the sign is kept 383 CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only 384 CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only 384 385 REAL(wp) , OPTIONAL , INTENT(in ) :: pval ! background value (used at closed boundaries) 385 386 !! … … 402 403 DO jk = 1, jpk 403 404 DO jj = nlcj+1, jpj ! added line(s) (inner only) 404 ptab(nldi :nlei , jj ,jk) = ptab(nldi:nlei, nlej,jk) 405 ptab(nldi :nlei , jj ,jk) = ptab(nldi:nlei, nlej,jk) 405 406 ptab(1 :nldi-1, jj ,jk) = ptab(nldi , nlej,jk) 406 407 ptab(nlei+1:nlci , jj ,jk) = ptab( nlei, nlej,jk) … … 413 414 END DO 414 415 ! 415 ELSE ! standard close or cyclic treatment 416 ELSE ! standard close or cyclic treatment 416 417 ! 417 418 ! ! East-West boundaries … … 432 433 ! 2. East and west directions exchange 433 434 ! ------------------------------------ 434 ! we play with the neigbours AND the row number because of the periodicity 435 ! we play with the neigbours AND the row number because of the periodicity 435 436 ! 436 437 SELECT CASE ( nbondi ) ! Read Dirichlet lateral conditions … … 441 442 t3we(:,jl,:,1) = ptab(iihom +jl,:,:) 442 443 END DO 443 END SELECT 444 END SELECT 444 445 ! 445 446 ! ! Migrations 446 447 imigr = jpreci * jpj * jpk 447 448 ! 448 SELECT CASE ( nbondi ) 449 SELECT CASE ( nbondi ) 449 450 CASE ( -1 ) 450 451 CALL mppsend( 2, t3we(1,1,1,1), imigr, noea, ml_req1 ) … … 472 473 ptab(iihom+jl,:,:) = t3ew(:,jl,:,2) 473 474 END DO 474 CASE ( 0 ) 475 CASE ( 0 ) 475 476 DO jl = 1, jpreci 476 477 ptab(jl ,:,:) = t3we(:,jl,:,2) … … 499 500 imigr = jprecj * jpi * jpk 500 501 ! 501 SELECT CASE ( nbondj ) 502 SELECT CASE ( nbondj ) 502 503 CASE ( -1 ) 503 504 CALL mppsend( 4, t3sn(1,1,1,1), imigr, nono, ml_req1 ) … … 511 512 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 512 513 IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 513 CASE ( 1 ) 514 CASE ( 1 ) 514 515 CALL mppsend( 3, t3ns(1,1,1,1), imigr, noso, ml_req1 ) 515 516 CALL mpprecv( 4, t3sn(1,1,1,2), imigr, noso ) … … 525 526 ptab(:,ijhom+jl,:) = t3ns(:,jl,:,2) 526 527 END DO 527 CASE ( 0 ) 528 CASE ( 0 ) 528 529 DO jl = 1, jprecj 529 530 ptab(:,jl ,:) = t3sn(:,jl,:,2) … … 555 556 !!---------------------------------------------------------------------- 556 557 !! *** routine mpp_lnk_2d *** 557 !! 558 !! 558 559 !! ** Purpose : Message passing manadgement for 2d array 559 560 !! 560 !! ** Method : Use mppsend and mpprecv function for passing mask 561 !! ** Method : Use mppsend and mpprecv function for passing mask 561 562 !! between processors following neighboring subdomains. 562 563 !! domain parameters … … 565 566 !! nbondi : mark for "east-west local boundary" 566 567 !! nbondj : mark for "north-south local boundary" 567 !! noea : number for local neighboring processors 568 !! noea : number for local neighboring processors 568 569 !! nowe : number for local neighboring processors 569 570 !! noso : number for local neighboring processors … … 576 577 REAL(wp) , INTENT(in ) :: psgn ! =-1 the sign change across the north fold boundary 577 578 ! ! = 1. , the sign is kept 578 CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only 579 CHARACTER(len=3), OPTIONAL , INTENT(in ) :: cd_mpp ! fill the overlap area only 579 580 REAL(wp) , OPTIONAL , INTENT(in ) :: pval ! background value (used at closed boundaries) 580 581 !! … … 597 598 ! WARNING pt2d is defined only between nld and nle 598 599 DO jj = nlcj+1, jpj ! added line(s) (inner only) 599 pt2d(nldi :nlei , jj ) = pt2d(nldi:nlei, nlej) 600 pt2d(nldi :nlei , jj ) = pt2d(nldi:nlei, nlej) 600 601 pt2d(1 :nldi-1, jj ) = pt2d(nldi , nlej) 601 602 pt2d(nlei+1:nlci , jj ) = pt2d( nlei, nlej) … … 607 608 END DO 608 609 ! 609 ELSE ! standard close or cyclic treatment 610 ELSE ! standard close or cyclic treatment 610 611 ! 611 612 ! ! East-West boundaries … … 626 627 ! 2. East and west directions exchange 627 628 ! ------------------------------------ 628 ! we play with the neigbours AND the row number because of the periodicity 629 ! we play with the neigbours AND the row number because of the periodicity 629 630 ! 630 631 SELECT CASE ( nbondi ) ! Read Dirichlet lateral conditions … … 724 725 pt2d(:,ijhom+jl) = t2ns(:,jl,2) 725 726 END DO 726 CASE ( 1 ) 727 CASE ( 1 ) 727 728 DO jl = 1, jprecj 728 729 pt2d(:,jl ) = t2sn(:,jl,2) … … 752 753 !! ** Purpose : Message passing manadgement for two 3D arrays 753 754 !! 754 !! ** Method : Use mppsend and mpprecv function for passing mask 755 !! ** Method : Use mppsend and mpprecv function for passing mask 755 756 !! between processors following neighboring subdomains. 756 757 !! domain parameters … … 759 760 !! nbondi : mark for "east-west local boundary" 760 761 !! nbondj : mark for "north-south local boundary" 761 !! noea : number for local neighboring processors 762 !! noea : number for local neighboring processors 762 763 !! nowe : number for local neighboring processors 763 764 !! noso : number for local neighboring processors … … 767 768 !! 768 769 !!---------------------------------------------------------------------- 769 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: ptab1 ! first and second 3D array on which 770 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: ptab1 ! first and second 3D array on which 770 771 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: ptab2 ! the boundary condition is applied 771 CHARACTER(len=1) , INTENT(in ) :: cd_type1 ! nature of ptab1 and ptab2 arrays 772 CHARACTER(len=1) , INTENT(in ) :: cd_type1 ! nature of ptab1 and ptab2 arrays 772 773 CHARACTER(len=1) , INTENT(in ) :: cd_type2 ! i.e. grid-points = T , U , V , F or W points 773 774 REAL(wp) , INTENT(in ) :: psgn ! =-1 the sign change across the north fold boundary … … 795 796 ENDIF 796 797 797 798 798 799 ! ! North-South boundaries 799 800 IF( .NOT. cd_type1 == 'F' ) ptab1(:, 1 :jprecj,:) = 0.e0 ! south except at F-point … … 805 806 ! 2. East and west directions exchange 806 807 ! ------------------------------------ 807 ! we play with the neigbours AND the row number because of the periodicity 808 ! we play with the neigbours AND the row number because of the periodicity 808 809 ! 809 810 SELECT CASE ( nbondi ) ! Read Dirichlet lateral conditions … … 821 822 imigr = jpreci * jpj * jpk *2 822 823 ! 823 SELECT CASE ( nbondi ) 824 SELECT CASE ( nbondi ) 824 825 CASE ( -1 ) 825 826 CALL mppsend( 2, t4we(1,1,1,1,1), imigr, noea, ml_req1 ) … … 848 849 ptab2(iihom+jl,:,:) = t4ew(:,jl,:,2,2) 849 850 END DO 850 CASE ( 0 ) 851 CASE ( 0 ) 851 852 DO jl = 1, jpreci 852 853 ptab1(jl ,:,:) = t4we(:,jl,:,1,2) … … 880 881 imigr = jprecj * jpi * jpk * 2 881 882 ! 882 SELECT CASE ( nbondj ) 883 SELECT CASE ( nbondj ) 883 884 CASE ( -1 ) 884 885 CALL mppsend( 4, t4sn(1,1,1,1,1), imigr, nono, ml_req1 ) … … 892 893 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 893 894 IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 894 CASE ( 1 ) 895 CASE ( 1 ) 895 896 CALL mppsend( 3, t4ns(1,1,1,1,1), imigr, noso, ml_req1 ) 896 897 CALL mpprecv( 4, t4sn(1,1,1,1,2), imigr, noso ) … … 907 908 ptab2(:,ijhom+jl,:) = t4ns(:,jl,:,2,2) 908 909 END DO 909 CASE ( 0 ) 910 CASE ( 0 ) 910 911 DO jl = 1, jprecj 911 912 ptab1(:,jl ,:) = t4sn(:,jl,:,1,2) … … 927 928 ! 928 929 SELECT CASE ( jpni ) 929 CASE ( 1 ) 930 CASE ( 1 ) 930 931 CALL lbc_nfd ( ptab1, cd_type1, psgn ) ! only for northern procs. 931 932 CALL lbc_nfd ( ptab2, cd_type2, psgn ) … … 933 934 CALL mpp_lbc_north( ptab1, cd_type1, psgn ) ! for all northern procs. 934 935 CALL mpp_lbc_north (ptab2, cd_type2, psgn) 935 END SELECT 936 END SELECT 936 937 ! 937 938 ENDIF … … 943 944 !!---------------------------------------------------------------------- 944 945 !! *** routine mpp_lnk_2d_e *** 945 !! 946 !! 946 947 !! ** Purpose : Message passing manadgement for 2d array (with halo) 947 948 !! 948 !! ** Method : Use mppsend and mpprecv function for passing mask 949 !! ** Method : Use mppsend and mpprecv function for passing mask 949 950 !! between processors following neighboring subdomains. 950 951 !! domain parameters … … 955 956 !! nbondi : mark for "east-west local boundary" 956 957 !! nbondj : mark for "north-south local boundary" 957 !! noea : number for local neighboring processors 958 !! noea : number for local neighboring processors 958 959 !! nowe : number for local neighboring processors 959 960 !! noso : number for local neighboring processors … … 984 985 IF( .NOT. cd_type == 'F' ) pt2d(:, 1-jpr2dj : jprecj ) = 0.e0 ! south except at F-point 985 986 pt2d(:,nlcj-jprecj+1:jpj+jpr2dj) = 0.e0 ! north 986 987 987 988 ! ! East-West boundaries 988 989 ! !* Cyclic east-west … … 1004 1005 CASE ( 1 ) ; CALL lbc_nfd ( pt2d(1:jpi,1:jpj+jpr2dj), cd_type, psgn, pr2dj=jpr2dj ) 1005 1006 CASE DEFAULT ; CALL mpp_lbc_north_e( pt2d , cd_type, psgn ) 1006 END SELECT 1007 END SELECT 1007 1008 ! 1008 1009 ENDIF … … 1010 1011 ! 2. East and west directions exchange 1011 1012 ! ------------------------------------ 1012 ! we play with the neigbours AND the row number because of the periodicity 1013 ! we play with the neigbours AND the row number because of the periodicity 1013 1014 ! 1014 1015 SELECT CASE ( nbondi ) ! Read Dirichlet lateral conditions … … 1096 1097 ! 1097 1098 ! ! Write Dirichlet lateral conditions 1098 ijhom = nlcj - jprecj 1099 ijhom = nlcj - jprecj 1099 1100 ! 1100 1101 SELECT CASE ( nbondj ) … … 1108 1109 pt2d(:,ijhom+jl ) = tr2ns(:,jl,2) 1109 1110 END DO 1110 CASE ( 1 ) 1111 CASE ( 1 ) 1111 1112 DO jl = 1, iprecj 1112 1113 pt2d(:,jl-jpr2dj) = tr2sn(:,jl,2) … … 1120 1121 !!---------------------------------------------------------------------- 1121 1122 !! *** routine mppsend *** 1122 !! 1123 !! 1123 1124 !! ** Purpose : Send messag passing array 1124 1125 !! … … 1156 1157 INTEGER , INTENT(in ) :: kbytes ! suze of the array pmess 1157 1158 INTEGER , INTENT(in ) :: ktyp ! Tag of the recevied message 1158 INTEGER, OPTIONAL, INTENT(in) :: ksource ! source process number 1159 INTEGER, OPTIONAL, INTENT(in) :: ksource ! source process number 1159 1160 !! 1160 1161 INTEGER :: istatus(mpi_status_size) … … 1164 1165 ! 1165 1166 1166 ! If a specific process number has been passed to the receive call, 1167 ! If a specific process number has been passed to the receive call, 1167 1168 ! use that one. Default is to use mpi_any_source 1168 1169 use_source=mpi_any_source … … 1179 1180 !!---------------------------------------------------------------------- 1180 1181 !! *** routine mppgather *** 1181 !! 1182 !! ** Purpose : Transfert between a local subdomain array and a work 1182 !! 1183 !! ** Purpose : Transfert between a local subdomain array and a work 1183 1184 !! array which is distributed following the vertical level. 1184 1185 !! … … 1193 1194 itaille = jpi * jpj 1194 1195 CALL mpi_gather( ptab, itaille, mpi_double_precision, pio, itaille , & 1195 & mpi_double_precision, kp , mpi_comm_opa, ierror ) 1196 & mpi_double_precision, kp , mpi_comm_opa, ierror ) 1196 1197 ! 1197 1198 END SUBROUTINE mppgather … … 1202 1203 !! *** routine mppscatter *** 1203 1204 !! 1204 !! ** Purpose : Transfert between awork array which is distributed 1205 !! ** Purpose : Transfert between awork array which is distributed 1205 1206 !! following the vertical level and the local subdomain array. 1206 1207 !! … … 1224 1225 !!---------------------------------------------------------------------- 1225 1226 !! *** routine mppmax_a_int *** 1226 !! 1227 !! 1227 1228 !! ** Purpose : Find maximum value in an integer layout array 1228 1229 !! … … 1230 1231 INTEGER , INTENT(in ) :: kdim ! size of array 1231 1232 INTEGER , INTENT(inout), DIMENSION(kdim) :: ktab ! input array 1232 INTEGER , INTENT(in ), OPTIONAL :: kcom ! 1233 INTEGER , INTENT(in ), OPTIONAL :: kcom ! 1233 1234 !! 1234 1235 INTEGER :: ierror, localcomm ! temporary integer … … 1255 1256 INTEGER, INTENT(inout) :: ktab ! ??? 1256 1257 INTEGER, INTENT(in ), OPTIONAL :: kcom ! ??? 1257 !! 1258 !! 1258 1259 INTEGER :: ierror, iwork, localcomm ! temporary integer 1259 1260 !!---------------------------------------------------------------------- 1260 1261 ! 1261 localcomm = mpi_comm_opa 1262 localcomm = mpi_comm_opa 1262 1263 IF( PRESENT(kcom) ) localcomm = kcom 1263 1264 ! … … 1272 1273 !!---------------------------------------------------------------------- 1273 1274 !! *** routine mppmin_a_int *** 1274 !! 1275 !! 1275 1276 !! ** Purpose : Find minimum value in an integer layout array 1276 1277 !! … … 1320 1321 !!---------------------------------------------------------------------- 1321 1322 !! *** routine mppsum_a_int *** 1322 !! 1323 !! 1323 1324 !! ** Purpose : Global integer sum, 1D array case 1324 1325 !! … … 1341 1342 !!---------------------------------------------------------------------- 1342 1343 !! *** routine mppsum_int *** 1343 !! 1344 !! 1344 1345 !! ** Purpose : Global integer sum 1345 1346 !! 1346 1347 !!---------------------------------------------------------------------- 1347 1348 INTEGER, INTENT(inout) :: ktab 1348 !! 1349 !! 1349 1350 INTEGER :: ierror, iwork 1350 1351 !!---------------------------------------------------------------------- … … 1360 1361 !!---------------------------------------------------------------------- 1361 1362 !! *** routine mppmax_a_real *** 1362 !! 1363 !! 1363 1364 !! ** Purpose : Maximum 1364 1365 !! … … 1384 1385 !!---------------------------------------------------------------------- 1385 1386 !! *** routine mppmax_real *** 1386 !! 1387 !! 1387 1388 !! ** Purpose : Maximum 1388 1389 !! … … 1395 1396 !!---------------------------------------------------------------------- 1396 1397 ! 1397 localcomm = mpi_comm_opa 1398 localcomm = mpi_comm_opa 1398 1399 IF( PRESENT(kcom) ) localcomm = kcom 1399 1400 ! … … 1407 1408 !!---------------------------------------------------------------------- 1408 1409 !! *** routine mppmin_a_real *** 1409 !! 1410 !! 1410 1411 !! ** Purpose : Minimum of REAL, array case 1411 1412 !! … … 1419 1420 !!----------------------------------------------------------------------- 1420 1421 ! 1421 localcomm = mpi_comm_opa 1422 localcomm = mpi_comm_opa 1422 1423 IF( PRESENT(kcom) ) localcomm = kcom 1423 1424 ! … … 1431 1432 !!---------------------------------------------------------------------- 1432 1433 !! *** routine mppmin_real *** 1433 !! 1434 !! 1434 1435 !! ** Purpose : minimum of REAL, scalar case 1435 1436 !! 1436 1437 !!----------------------------------------------------------------------- 1437 REAL(wp), INTENT(inout) :: ptab ! 1438 REAL(wp), INTENT(inout) :: ptab ! 1438 1439 INTEGER , INTENT(in ), OPTIONAL :: kcom 1439 1440 !! … … 1443 1444 !!----------------------------------------------------------------------- 1444 1445 ! 1445 localcomm = mpi_comm_opa 1446 localcomm = mpi_comm_opa 1446 1447 IF( PRESENT(kcom) ) localcomm = kcom 1447 1448 ! … … 1455 1456 !!---------------------------------------------------------------------- 1456 1457 !! *** routine mppsum_a_real *** 1457 !! 1458 !! 1458 1459 !! ** Purpose : global sum, REAL ARRAY argument case 1459 1460 !! … … 1464 1465 !! 1465 1466 INTEGER :: ierror ! temporary integer 1466 INTEGER :: localcomm 1467 REAL(wp), DIMENSION(kdim) :: zwork ! temporary workspace 1467 INTEGER :: localcomm 1468 REAL(wp), DIMENSION(kdim) :: zwork ! temporary workspace 1468 1469 !!----------------------------------------------------------------------- 1469 1470 ! 1470 localcomm = mpi_comm_opa 1471 localcomm = mpi_comm_opa 1471 1472 IF( PRESENT(kcom) ) localcomm = kcom 1472 1473 ! … … 1480 1481 !!---------------------------------------------------------------------- 1481 1482 !! *** routine mppsum_real *** 1482 !! 1483 !! 1483 1484 !! ** Purpose : global sum, SCALAR argument case 1484 1485 !! … … 1487 1488 INTEGER , INTENT(in ), OPTIONAL :: kcom 1488 1489 !! 1489 INTEGER :: ierror, localcomm 1490 INTEGER :: ierror, localcomm 1490 1491 REAL(wp) :: zwork 1491 1492 !!----------------------------------------------------------------------- 1492 1493 ! 1493 localcomm = mpi_comm_opa 1494 localcomm = mpi_comm_opa 1494 1495 IF( PRESENT(kcom) ) localcomm = kcom 1495 1496 ! … … 1524 1525 1525 1526 END SUBROUTINE mppsum_realdd 1526 1527 1527 1528 1528 1529 SUBROUTINE mppsum_a_realdd( ytab, kdim, kcom ) 1529 1530 !!---------------------------------------------------------------------- … … 1551 1552 1552 1553 END SUBROUTINE mppsum_a_realdd 1553 1554 1554 1555 SUBROUTINE mpp_minloc2d( ptab, pmask, pmin, ki,kj ) 1555 1556 !!------------------------------------------------------------------------ … … 1646 1647 REAL(wp) , INTENT( out) :: pmax ! Global maximum of ptab 1647 1648 INTEGER , INTENT( out) :: ki, kj ! index of maximum in global frame 1648 !! 1649 !! 1649 1650 INTEGER :: ierror 1650 1651 INTEGER, DIMENSION (2) :: ilocs … … 1685 1686 REAL(wp) , INTENT( out) :: pmax ! Global maximum of ptab 1686 1687 INTEGER , INTENT( out) :: ki, kj, kk ! index of maximum in global frame 1687 !! 1688 !! 1688 1689 REAL(wp) :: zmax ! local maximum 1689 1690 REAL(wp), DIMENSION(2,1) :: zain, zaout … … 1715 1716 !!---------------------------------------------------------------------- 1716 1717 !! *** routine mppsync *** 1717 !! 1718 !! 1718 1719 !! ** Purpose : Massively parallel processors, synchroneous 1719 1720 !! … … 1730 1731 !!---------------------------------------------------------------------- 1731 1732 !! *** routine mppstop *** 1732 !! 1733 !! 1733 1734 !! ** purpose : Stop massively parallel processors method 1734 1735 !! … … 1746 1747 !!---------------------------------------------------------------------- 1747 1748 !! *** routine mppobc *** 1748 !! 1749 !! 1749 1750 !! ** Purpose : Message passing manadgement for open boundary 1750 1751 !! conditions array … … 1757 1758 !! nbondi : mark for "east-west local boundary" 1758 1759 !! nbondj : mark for "north-south local boundary" 1759 !! noea : number for local neighboring processors 1760 !! noea : number for local neighboring processors 1760 1761 !! nowe : number for local neighboring processors 1761 1762 !! noso : number for local neighboring processors … … 1806 1807 CALL mppstop 1807 1808 ENDIF 1808 1809 1809 1810 ! Communication level by level 1810 1811 ! ---------------------------- … … 1921 1922 DO jj = ijpt0, ijpt1 ! north/south boundaries 1922 1923 DO ji = iipt0,ilpt1 1923 ptab(ji,jk) = ztab(ji,jj) 1924 ptab(ji,jk) = ztab(ji,jj) 1924 1925 END DO 1925 1926 END DO … … 1927 1928 DO jj = ijpt0, ilpt1 ! east/west boundaries 1928 1929 DO ji = iipt0,iipt1 1929 ptab(jj,jk) = ztab(ji,jj) 1930 ptab(jj,jk) = ztab(ji,jj) 1930 1931 END DO 1931 1932 END DO … … 1937 1938 ! 1938 1939 END SUBROUTINE mppobc 1939 1940 1940 1941 1941 1942 SUBROUTINE mpp_comm_free( kcom ) … … 1996 1997 kice = 0 1997 1998 DO jjproc = 1, jpnij 1998 IF( jjproc == narea .AND. pindic .GT. 0 ) kice(jjproc) = 1 1999 IF( jjproc == narea .AND. pindic .GT. 0 ) kice(jjproc) = 1 1999 2000 END DO 2000 2001 ! 2001 2002 zwork = 0 2002 2003 CALL MPI_ALLREDUCE( kice, zwork, jpnij, mpi_integer, mpi_sum, mpi_comm_opa, ierr ) 2003 ndim_rank_ice = SUM( zwork ) 2004 ndim_rank_ice = SUM( zwork ) 2004 2005 2005 2006 ! Allocate the right size to nrank_north … … 2007 2008 ALLOCATE( nrank_ice(ndim_rank_ice) ) 2008 2009 ! 2009 ii = 0 2010 ii = 0 2010 2011 nrank_ice = 0 2011 2012 DO jjproc = 1, jpnij 2012 2013 IF( zwork(jjproc) == 1) THEN 2013 2014 ii = ii + 1 2014 nrank_ice(ii) = jjproc -1 2015 nrank_ice(ii) = jjproc -1 2015 2016 ENDIF 2016 2017 END DO … … 2094 2095 IF (ALLOCATED (nrank_znl)) DEALLOCATE(nrank_znl) 2095 2096 ALLOCATE(nrank_znl(ndim_rank_znl)) 2096 ii = 0 2097 ii = 0 2097 2098 nrank_znl (:) = 0 2098 2099 DO jproc=1,jpnij 2099 2100 IF ( kwork(jproc) == njmpp) THEN 2100 2101 ii = ii + 1 2101 nrank_znl(ii) = jproc -1 2102 nrank_znl(ii) = jproc -1 2102 2103 ENDIF 2103 2104 END DO … … 2123 2124 2124 2125 ! Determines if processor if the first (starting from i=1) on the row 2125 IF ( jpni == 1 ) THEN 2126 IF ( jpni == 1 ) THEN 2126 2127 l_znl_root = .TRUE. 2127 2128 ELSE … … 2141 2142 !! *** routine mpp_ini_north *** 2142 2143 !! 2143 !! ** Purpose : Initialize special communicator for north folding 2144 !! ** Purpose : Initialize special communicator for north folding 2144 2145 !! condition together with global variables needed in the mpp folding 2145 2146 !! … … 2202 2203 !! *** routine mpp_lbc_north_3d *** 2203 2204 !! 2204 !! ** Purpose : Ensure proper north fold horizontal bondary condition 2205 !! ** Purpose : Ensure proper north fold horizontal bondary condition 2205 2206 !! in mpp configuration in case of jpn1 > 1 2206 2207 !! 2207 2208 !! ** Method : North fold condition and mpp with more than one proc 2208 !! in i-direction require a specific treatment. We gather 2209 !! in i-direction require a specific treatment. We gather 2209 2210 !! the 4 northern lines of the global domain on 1 processor 2210 2211 !! and apply lbc north-fold on this sub array. Then we … … 2215 2216 CHARACTER(len=1) , INTENT(in ) :: cd_type ! nature of pt3d grid-points 2216 2217 ! ! = T , U , V , F or W gridpoints 2217 REAL(wp) , INTENT(in ) :: psgn ! = -1. the sign change across the north fold 2218 REAL(wp) , INTENT(in ) :: psgn ! = -1. the sign change across the north fold 2218 2219 !! ! = 1. , the sign is kept 2219 2220 INTEGER :: ji, jj, jr … … 2224 2225 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat ! for mpi_isend when avoiding mpi_allgather 2225 2226 !!---------------------------------------------------------------------- 2226 ! 2227 ! 2227 2228 ijpj = 4 2228 2229 ityp = -1 … … 2239 2240 IF ( l_north_nogather ) THEN 2240 2241 ! 2241 ! Avoid the use of mpi_allgather by exchanging only with the processes already identified 2242 ! Avoid the use of mpi_allgather by exchanging only with the processes already identified 2242 2243 ! (in nemo_northcomms) as being involved in this process' northern boundary exchange 2243 2244 ! … … 2264 2265 ityp = 5 2265 2266 CASE DEFAULT 2266 ityp = -1 ! Set a default value for unsupported types which 2267 ityp = -1 ! Set a default value for unsupported types which 2267 2268 ! will cause a fallback to the mpi_allgather method 2268 2269 END SELECT … … 2313 2314 ! The ztab array has been either: 2314 2315 ! a. Fully populated by the mpi_allgather operation or 2315 ! b. Had the active points for this domain and northern neighbours populated 2316 ! b. Had the active points for this domain and northern neighbours populated 2316 2317 ! by peer to peer exchanges 2317 ! Either way the array may be folded by lbc_nfd and the result for the span of 2318 ! Either way the array may be folded by lbc_nfd and the result for the span of 2318 2319 ! this domain will be identical. 2319 2320 ! … … 2334 2335 !! *** routine mpp_lbc_north_2d *** 2335 2336 !! 2336 !! ** Purpose : Ensure proper north fold horizontal bondary condition 2337 !! ** Purpose : Ensure proper north fold horizontal bondary condition 2337 2338 !! in mpp configuration in case of jpn1 > 1 (for 2d array ) 2338 2339 !! 2339 2340 !! ** Method : North fold condition and mpp with more than one proc 2340 !! in i-direction require a specific treatment. We gather 2341 !! in i-direction require a specific treatment. We gather 2341 2342 !! the 4 northern lines of the global domain on 1 processor 2342 2343 !! and apply lbc north-fold on this sub array. Then we … … 2347 2348 CHARACTER(len=1) , INTENT(in ) :: cd_type ! nature of pt3d grid-points 2348 2349 ! ! = T , U , V , F or W gridpoints 2349 REAL(wp) , INTENT(in ) :: psgn ! = -1. the sign change across the north fold 2350 REAL(wp) , INTENT(in ) :: psgn ! = -1. the sign change across the north fold 2350 2351 !! ! = 1. , the sign is kept 2351 2352 INTEGER :: ji, jj, jr … … 2371 2372 IF ( l_north_nogather ) THEN 2372 2373 ! 2373 ! Avoid the use of mpi_allgather by exchanging only with the processes already identified 2374 ! Avoid the use of mpi_allgather by exchanging only with the processes already identified 2374 2375 ! (in nemo_northcomms) as being involved in this process' northern boundary exchange 2375 2376 ! … … 2396 2397 ityp = 5 2397 2398 CASE DEFAULT 2398 ityp = -1 ! Set a default value for unsupported types which 2399 ityp = -1 ! Set a default value for unsupported types which 2399 2400 ! will cause a fallback to the mpi_allgather method 2400 2401 END SELECT … … 2446 2447 ! The ztab array has been either: 2447 2448 ! a. Fully populated by the mpi_allgather operation or 2448 ! b. Had the active points for this domain and northern neighbours populated 2449 ! b. Had the active points for this domain and northern neighbours populated 2449 2450 ! by peer to peer exchanges 2450 ! Either way the array may be folded by lbc_nfd and the result for the span of 2451 ! Either way the array may be folded by lbc_nfd and the result for the span of 2451 2452 ! this domain will be identical. 2452 2453 ! … … 2468 2469 !! *** routine mpp_lbc_north_2d *** 2469 2470 !! 2470 !! ** Purpose : Ensure proper north fold horizontal bondary condition 2471 !! in mpp configuration in case of jpn1 > 1 and for 2d 2471 !! ** Purpose : Ensure proper north fold horizontal bondary condition 2472 !! in mpp configuration in case of jpn1 > 1 and for 2d 2472 2473 !! array with outer extra halo 2473 2474 !! 2474 2475 !! ** Method : North fold condition and mpp with more than one proc 2475 !! in i-direction require a specific treatment. We gather 2476 !! the 4+2*jpr2dj northern lines of the global domain on 1 2477 !! processor and apply lbc north-fold on this sub array. 2476 !! in i-direction require a specific treatment. We gather 2477 !! the 4+2*jpr2dj northern lines of the global domain on 1 2478 !! processor and apply lbc north-fold on this sub array. 2478 2479 !! Then we scatter the north fold array back to the processors. 2479 2480 !! … … 2482 2483 CHARACTER(len=1) , INTENT(in ) :: cd_type ! nature of pt3d grid-points 2483 2484 ! ! = T , U , V , F or W -points 2484 REAL(wp) , INTENT(in ) :: psgn ! = -1. the sign change across the 2485 REAL(wp) , INTENT(in ) :: psgn ! = -1. the sign change across the 2485 2486 !! ! north fold, = 1. otherwise 2486 2487 INTEGER :: ji, jj, jr … … 2525 2526 !! Scatter back to pt2d 2526 2527 DO jj = nlcj - ijpj + 1 , nlcj +jpr2dj 2527 ij = ij +1 2528 ij = ij +1 2528 2529 DO ji= 1, nlci 2529 2530 pt2d(ji,jj) = ztab_e(ji+nimpp-1,ij) … … 2542 2543 !! ** Method :: define buffer size in namelist, if 0 no buffer attachment 2543 2544 !! but classical mpi_init 2544 !! 2545 !! History :: 01/11 :: IDRIS initial version for IBM only 2545 !! 2546 !! History :: 01/11 :: IDRIS initial version for IBM only 2546 2547 !! 08/04 :: R. Benshila, generalisation 2547 2548 !!--------------------------------------------------------------------- 2548 CHARACTER(len=*),DIMENSION(:), INTENT( out) :: ldtxt 2549 CHARACTER(len=*),DIMENSION(:), INTENT( out) :: ldtxt 2549 2550 INTEGER , INTENT(inout) :: ksft 2550 2551 INTEGER , INTENT( out) :: code … … 2555 2556 CALL mpi_initialized( mpi_was_called, code ) ! MPI initialization 2556 2557 IF ( code /= MPI_SUCCESS ) THEN 2557 DO ji = 1, SIZE(ldtxt) 2558 DO ji = 1, SIZE(ldtxt) 2558 2559 IF( TRIM(ldtxt(ji)) /= '' ) WRITE(*,*) ldtxt(ji) ! control print of mynode 2559 END DO 2560 END DO 2560 2561 WRITE(*, cform_err) 2561 2562 WRITE(*, *) ' lib_mpp: Error in routine mpi_initialized' … … 2567 2568 CALL mpi_comm_dup( mpi_comm_world, mpi_comm_opa, code ) 2568 2569 IF ( code /= MPI_SUCCESS ) THEN 2569 DO ji = 1, SIZE(ldtxt) 2570 DO ji = 1, SIZE(ldtxt) 2570 2571 IF( TRIM(ldtxt(ji)) /= '' ) WRITE(*,*) ldtxt(ji) ! control print of mynode 2571 2572 END DO … … 2580 2581 ! Buffer allocation and attachment 2581 2582 ALLOCATE( tampon(nn_buffer), stat = ierr ) 2582 IF( ierr /= 0 ) THEN 2583 DO ji = 1, SIZE(ldtxt) 2583 IF( ierr /= 0 ) THEN 2584 DO ji = 1, SIZE(ldtxt) 2584 2585 IF( TRIM(ldtxt(ji)) /= '' ) WRITE(*,*) ldtxt(ji) ! control print of mynode 2585 2586 END DO … … 2660 2661 FUNCTION mynode( ldtxt, kumnam, kstop, localComm ) RESULT (function_value) 2661 2662 INTEGER, OPTIONAL , INTENT(in ) :: localComm 2662 CHARACTER(len=*),DIMENSION(:) :: ldtxt 2663 CHARACTER(len=*),DIMENSION(:) :: ldtxt 2663 2664 INTEGER :: kumnam, kstop 2664 2665 IF( PRESENT( localComm ) .OR. .NOT.PRESENT( localComm ) ) function_value = 0 … … 2672 2673 REAL , DIMENSION(:) :: parr 2673 2674 INTEGER :: kdim 2674 INTEGER, OPTIONAL :: kcom 2675 INTEGER, OPTIONAL :: kcom 2675 2676 WRITE(*,*) 'mpp_sum_as: You should not have seen this print! error?', kdim, parr(1), kcom 2676 2677 END SUBROUTINE mpp_sum_as … … 2679 2680 REAL , DIMENSION(:,:) :: parr 2680 2681 INTEGER :: kdim 2681 INTEGER, OPTIONAL :: kcom 2682 INTEGER, OPTIONAL :: kcom 2682 2683 WRITE(*,*) 'mpp_sum_a2s: You should not have seen this print! error?', kdim, parr(1,1), kcom 2683 2684 END SUBROUTINE mpp_sum_a2s … … 2686 2687 INTEGER, DIMENSION(:) :: karr 2687 2688 INTEGER :: kdim 2688 INTEGER, OPTIONAL :: kcom 2689 INTEGER, OPTIONAL :: kcom 2689 2690 WRITE(*,*) 'mpp_sum_ai: You should not have seen this print! error?', kdim, karr(1), kcom 2690 2691 END SUBROUTINE mpp_sum_ai … … 2692 2693 SUBROUTINE mpp_sum_s( psca, kcom ) ! Dummy routine 2693 2694 REAL :: psca 2694 INTEGER, OPTIONAL :: kcom 2695 INTEGER, OPTIONAL :: kcom 2695 2696 WRITE(*,*) 'mpp_sum_s: You should not have seen this print! error?', psca, kcom 2696 2697 END SUBROUTINE mpp_sum_s … … 2698 2699 SUBROUTINE mpp_sum_i( kint, kcom ) ! Dummy routine 2699 2700 integer :: kint 2700 INTEGER, OPTIONAL :: kcom 2701 INTEGER, OPTIONAL :: kcom 2701 2702 WRITE(*,*) 'mpp_sum_i: You should not have seen this print! error?', kint, kcom 2702 2703 END SUBROUTINE mpp_sum_i … … 2707 2708 WRITE(*,*) 'mppsum_realdd: You should not have seen this print! error?', ytab 2708 2709 END SUBROUTINE mppsum_realdd 2709 2710 2710 2711 SUBROUTINE mppsum_a_realdd( ytab, kdim, kcom ) 2711 2712 INTEGER , INTENT( in ) :: kdim ! size of ytab … … 2718 2719 REAL , DIMENSION(:) :: parr 2719 2720 INTEGER :: kdim 2720 INTEGER, OPTIONAL :: kcom 2721 INTEGER, OPTIONAL :: kcom 2721 2722 WRITE(*,*) 'mppmax_a_real: You should not have seen this print! error?', kdim, parr(1), kcom 2722 2723 END SUBROUTINE mppmax_a_real … … 2724 2725 SUBROUTINE mppmax_real( psca, kcom ) 2725 2726 REAL :: psca 2726 INTEGER, OPTIONAL :: kcom 2727 INTEGER, OPTIONAL :: kcom 2727 2728 WRITE(*,*) 'mppmax_real: You should not have seen this print! error?', psca, kcom 2728 2729 END SUBROUTINE mppmax_real … … 2731 2732 REAL , DIMENSION(:) :: parr 2732 2733 INTEGER :: kdim 2733 INTEGER, OPTIONAL :: kcom 2734 INTEGER, OPTIONAL :: kcom 2734 2735 WRITE(*,*) 'mppmin_a_real: You should not have seen this print! error?', kdim, parr(1), kcom 2735 2736 END SUBROUTINE mppmin_a_real … … 2737 2738 SUBROUTINE mppmin_real( psca, kcom ) 2738 2739 REAL :: psca 2739 INTEGER, OPTIONAL :: kcom 2740 INTEGER, OPTIONAL :: kcom 2740 2741 WRITE(*,*) 'mppmin_real: You should not have seen this print! error?', psca, kcom 2741 2742 END SUBROUTINE mppmin_real … … 2744 2745 INTEGER, DIMENSION(:) :: karr 2745 2746 INTEGER :: kdim 2746 INTEGER, OPTIONAL :: kcom 2747 INTEGER, OPTIONAL :: kcom 2747 2748 WRITE(*,*) 'mppmax_a_int: You should not have seen this print! error?', kdim, karr(1), kcom 2748 2749 END SUBROUTINE mppmax_a_int … … 2750 2751 SUBROUTINE mppmax_int( kint, kcom) 2751 2752 INTEGER :: kint 2752 INTEGER, OPTIONAL :: kcom 2753 INTEGER, OPTIONAL :: kcom 2753 2754 WRITE(*,*) 'mppmax_int: You should not have seen this print! error?', kint, kcom 2754 2755 END SUBROUTINE mppmax_int … … 2757 2758 INTEGER, DIMENSION(:) :: karr 2758 2759 INTEGER :: kdim 2759 INTEGER, OPTIONAL :: kcom 2760 INTEGER, OPTIONAL :: kcom 2760 2761 WRITE(*,*) 'mppmin_a_int: You should not have seen this print! error?', kdim, karr(1), kcom 2761 2762 END SUBROUTINE mppmin_a_int … … 2763 2764 SUBROUTINE mppmin_int( kint, kcom ) 2764 2765 INTEGER :: kint 2765 INTEGER, OPTIONAL :: kcom 2766 INTEGER, OPTIONAL :: kcom 2766 2767 WRITE(*,*) 'mppmin_int: You should not have seen this print! error?', kint, kcom 2767 2768 END SUBROUTINE mppmin_int … … 2850 2851 !! *** ROUTINE stop_opa *** 2851 2852 !! 2852 !! ** Purpose : print in ocean.outpput file a error message and 2853 !! ** Purpose : print in ocean.outpput file a error message and 2853 2854 !! increment the error number (nstop) by one. 2854 2855 !!---------------------------------------------------------------------- … … 2857 2858 !!---------------------------------------------------------------------- 2858 2859 ! 2859 nstop = nstop + 1 2860 nstop = nstop + 1 2860 2861 IF(lwp) THEN 2861 2862 WRITE(numout,cform_err) … … 2889 2890 !! *** ROUTINE stop_warn *** 2890 2891 !! 2891 !! ** Purpose : print in ocean.outpput file a error message and 2892 !! ** Purpose : print in ocean.outpput file a error message and 2892 2893 !! increment the warning number (nwarn) by one. 2893 2894 !!---------------------------------------------------------------------- … … 2895 2896 CHARACTER(len=*), INTENT(in), OPTIONAL :: cd6, cd7, cd8, cd9, cd10 2896 2897 !!---------------------------------------------------------------------- 2897 ! 2898 nwarn = nwarn + 1 2898 ! 2899 nwarn = nwarn + 1 2899 2900 IF(lwp) THEN 2900 2901 WRITE(numout,cform_war) … … 2982 2983 STOP 'ctl_opn bad opening' 2983 2984 ENDIF 2984 2985 2985 2986 END SUBROUTINE ctl_opn 2986 2987 … … 2992 2993 !! ** Purpose : return the index of an unused logical unit 2993 2994 !!---------------------------------------------------------------------- 2994 LOGICAL :: llopn 2995 LOGICAL :: llopn 2995 2996 !!---------------------------------------------------------------------- 2996 2997 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
r3421 r3598 6 6 !! History : OPA ! 2000-11 (R. Hordoir, E. Durand) NetCDF FORMAT 7 7 !! NEMO 1.0 ! 2002-09 (G. Madec) F90: Free form and module 8 !! 3.0 ! 2006-07 (G. Madec) Surface module 8 !! 3.0 ! 2006-07 (G. Madec) Surface module 9 9 !! 3.2 ! 2009-04 (B. Lemaire) Introduce iom_put 10 10 !! 3.3 ! 2010-10 (R. Furner, G. Madec) runoff distributed over ocean levels … … 32 32 PUBLIC sbc_rnf_div ! routine called in sshwzv module 33 33 PUBLIC sbc_rnf_alloc ! routine call in sbcmod module 34 34 PUBLIC sbc_rnf_init ! (PUBLIC for TAM) 35 35 ! !!* namsbc_rnf namelist * 36 36 CHARACTER(len=100), PUBLIC :: cn_dir = './' !: Root directory for location of ssr files 37 37 LOGICAL , PUBLIC :: ln_rnf_depth = .false. !: depth river runoffs attribute specified in a file 38 LOGICAL , PUBLIC :: ln_rnf_tem = .false. !: temperature river runoffs attribute specified in a file 39 LOGICAL , PUBLIC :: ln_rnf_sal = .false. !: salinity river runoffs attribute specified in a file 38 LOGICAL , PUBLIC :: ln_rnf_tem = .false. !: temperature river runoffs attribute specified in a file 39 LOGICAL , PUBLIC :: ln_rnf_sal = .false. !: salinity river runoffs attribute specified in a file 40 40 LOGICAL , PUBLIC :: ln_rnf_emp = .false. !: runoffs into a file to be read or already into precipitation 41 41 TYPE(FLD_N) , PUBLIC :: sn_rnf !: information about the runoff file to be read 42 42 TYPE(FLD_N) , PUBLIC :: sn_cnf !: information about the runoff mouth file to be read 43 TYPE(FLD_N) :: sn_s_rnf !: information about the salinities of runoff file to be read 44 TYPE(FLD_N) :: sn_t_rnf !: information about the temperatures of runoff file to be read 43 TYPE(FLD_N) :: sn_s_rnf !: information about the salinities of runoff file to be read 44 TYPE(FLD_N) :: sn_t_rnf !: information about the temperatures of runoff file to be read 45 45 TYPE(FLD_N) :: sn_dep_rnf !: information about the depth which river inflow affects 46 46 LOGICAL , PUBLIC :: ln_rnf_mouth = .false. !: specific treatment in mouths vicinity … … 55 55 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nk_rnf !: depth of runoff in model levels 56 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rnf_tsc_b, rnf_tsc !: before and now T & S runoff contents [K.m/s & PSU.m/s] 57 58 REAL(wp) :: r1_rau0 ! = 1 / rau0 59 60 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnf ! structure: river runoff (file information, fields read)61 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_s_rnf ! structure: river runoff salinity (file information, fields read)62 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_t_rnf ! structure: river runoff temperature (file information, fields read)63 64 !! * Substitutions 65 # include "domzgr_substitute.h90" 57 58 REAL(wp) :: r1_rau0 ! = 1 / rau0 59 60 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf_rnf ! structure: river runoff (file information, fields read) (PUBLIC for TAM) 61 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf_s_rnf ! structure: river runoff salinity (file information, fields read) (PUBLIC for TAM) 62 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf_t_rnf ! structure: river runoff temperature (file information, fields read) (PUBLIC for TAM) 63 64 !! * Substitutions 65 # include "domzgr_substitute.h90" 66 66 !!---------------------------------------------------------------------- 67 67 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 86 86 !!---------------------------------------------------------------------- 87 87 !! *** ROUTINE sbc_rnf *** 88 !! 88 !! 89 89 !! ** Purpose : Introduce a climatological run off forcing 90 90 !! 91 !! ** Method : Set each river mouth with a monthly climatology 91 !! ** Method : Set each river mouth with a monthly climatology 92 92 !! provided from different data. 93 93 !! CAUTION : upward water flux, runoff forced to be < 0 … … 99 99 INTEGER :: ji, jj ! dummy loop indices 100 100 !!---------------------------------------------------------------------- 101 ! 101 ! 102 102 IF( kt == nit000 ) CALL sbc_rnf_init ! Read namelist and allocate structures 103 103 … … 114 114 ! !-------------------! 115 115 ! 116 CALL fld_read ( kt, nn_fsbc, sf_rnf ) ! Read Runoffs data and provide it at kt 116 CALL fld_read ( kt, nn_fsbc, sf_rnf ) ! Read Runoffs data and provide it at kt 117 117 IF( ln_rnf_tem ) CALL fld_read ( kt, nn_fsbc, sf_t_rnf ) ! idem for runoffs temperature if required 118 118 IF( ln_rnf_sal ) CALL fld_read ( kt, nn_fsbc, sf_s_rnf ) ! idem for runoffs salinity if required … … 127 127 ! 128 128 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 129 rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) 129 rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) 130 130 ! 131 131 r1_rau0 = 1._wp / rau0 … … 133 133 IF( ln_rnf_tem ) THEN ! use runoffs temperature data 134 134 rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 135 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999 ) ! if missing data value use SST as runoffs temperature 135 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999 ) ! if missing data value use SST as runoffs temperature 136 136 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 137 137 END WHERE 138 138 ELSE ! use SST as runoffs temperature 139 139 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 140 ENDIF 141 ! ! use runoffs salinity data 140 ENDIF 141 ! ! use runoffs salinity data 142 142 IF( ln_rnf_sal ) rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 143 143 ! ! else use S=0 for runoffs (done one for all in the init) 144 144 ! 145 145 IF( ln_rnf_tem .OR. ln_rnf_sal ) THEN ! runoffs as outflow: use ocean SST and SSS 146 WHERE( rnf(:,:) < 0._wp ) ! example baltic model when flow is out of domain 146 WHERE( rnf(:,:) < 0._wp ) ! example baltic model when flow is out of domain 147 147 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 148 148 rnf_tsc(:,:,jp_sal) = sss_m(:,:) * rnf(:,:) * r1_rau0 … … 158 158 ! ! ---------------------------------------- ! 159 159 IF( ln_rstart .AND. & !* Restart: read in restart file 160 & iom_varid( numror, 'rnf_b', ldstop = .FALSE. ) > 0 ) THEN 160 & iom_varid( numror, 'rnf_b', ldstop = .FALSE. ) > 0 ) THEN 161 161 IF(lwp) WRITE(numout,*) ' nit000-1 runoff forcing fields red in the restart file' 162 162 CALL iom_get( numror, jpdom_autoglo, 'rnf_b', rnf_b ) ! before runoff … … 165 165 ELSE !* no restart: set from nit000 values 166 166 IF(lwp) WRITE(numout,*) ' nit000-1 runoff forcing fields set to nit000' 167 rnf_b (:,: ) = rnf (:,: ) 168 rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) 167 rnf_b (:,: ) = rnf (:,: ) 168 rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) 169 169 ENDIF 170 170 ENDIF … … 187 187 !!---------------------------------------------------------------------- 188 188 !! *** ROUTINE sbc_rnf *** 189 !! 189 !! 190 190 !! ** Purpose : update the horizontal divergence with the runoff inflow 191 191 !! 192 !! ** Method : 193 !! CAUTION : rnf is positive (inflow) decreasing the 192 !! ** Method : 193 !! CAUTION : rnf is positive (inflow) decreasing the 194 194 !! divergence and expressed in m/s 195 195 !! … … 207 207 r1_rau0 = 1._wp / rau0 208 208 IF( ln_rnf_depth ) THEN !== runoff distributed over several levels ==! 209 IF( lk_vvl ) THEN ! variable volume case 209 IF( lk_vvl ) THEN ! variable volume case 210 210 DO jj = 1, jpj ! update the depth over which runoffs are distributed 211 211 DO ji = 1, jpi 212 h_rnf(ji,jj) = 0._wp 212 h_rnf(ji,jj) = 0._wp 213 213 DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres 214 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) ! to the bottom of the relevant grid box 215 END DO 214 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) ! to the bottom of the relevant grid box 215 END DO 216 216 ! ! apply the runoff input flow 217 217 DO jk = 1, nk_rnf(ji,jj) … … 249 249 !! ** Action : - read parameters 250 250 !!---------------------------------------------------------------------- 251 CHARACTER(len=32) :: rn_dep_file ! runoff file name 251 CHARACTER(len=32) :: rn_dep_file ! runoff file name 252 252 INTEGER :: ji, jj, jk ! dummy loop indices 253 253 INTEGER :: ierror, inum ! temporary integer 254 !! 254 !! 255 255 NAMELIST/namsbc_rnf/ cn_dir, ln_rnf_emp, ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, & 256 & sn_rnf, sn_cnf , sn_s_rnf , sn_t_rnf , sn_dep_rnf, & 257 & ln_rnf_mouth , rn_hrnf , rn_avt_rnf, rn_rfact 256 & sn_rnf, sn_cnf , sn_s_rnf , sn_t_rnf , sn_dep_rnf, & 257 & ln_rnf_mouth , rn_hrnf , rn_avt_rnf, rn_rfact 258 258 !!---------------------------------------------------------------------- 259 259 … … 267 267 sn_cnf = FLD_N( 'runoffs', 0 , 'sorunoff' , .FALSE. , .true. , 'yearly' , '' , '' ) 268 268 269 sn_s_rnf = FLD_N( 'runoffs', 24. , 'rosaline' , .TRUE. , .true. , 'yearly' , '' , '' ) 270 sn_t_rnf = FLD_N( 'runoffs', 24. , 'rotemper' , .TRUE. , .true. , 'yearly' , '' , '' ) 271 sn_dep_rnf = FLD_N( 'runoffs', 0. , 'rodepth' , .FALSE. , .true. , 'yearly' , '' , '' ) 269 sn_s_rnf = FLD_N( 'runoffs', 24. , 'rosaline' , .TRUE. , .true. , 'yearly' , '' , '' ) 270 sn_t_rnf = FLD_N( 'runoffs', 24. , 'rotemper' , .TRUE. , .true. , 'yearly' , '' , '' ) 271 sn_dep_rnf = FLD_N( 'runoffs', 0. , 'rodepth' , .FALSE. , .true. , 'yearly' , '' , '' ) 272 272 ! 273 273 REWIND ( numnam ) ! Read Namelist namsbc_rnf … … 284 284 WRITE(numout,*) ' river mouth additional Kz rn_avt_rnf = ', rn_avt_rnf 285 285 WRITE(numout,*) ' depth of river mouth additional mixing rn_hrnf = ', rn_hrnf 286 WRITE(numout,*) ' multiplicative factor for runoff rn_rfact = ', rn_rfact 286 WRITE(numout,*) ' multiplicative factor for runoff rn_rfact = ', rn_rfact 287 287 ENDIF 288 288 … … 297 297 IF(lwp) WRITE(numout,*) ' runoffs directly provided in the precipitations' 298 298 IF( ln_rnf_depth .OR. ln_rnf_tem .OR. ln_rnf_sal ) THEN 299 CALL ctl_warn( 'runoffs already included in precipitations, so runoff (T,S, depth) attributes will not be used' ) 299 CALL ctl_warn( 'runoffs already included in precipitations, so runoff (T,S, depth) attributes will not be used' ) 300 300 ln_rnf_depth = .FALSE. ; ln_rnf_tem = .FALSE. ; ln_rnf_sal = .FALSE. 301 301 ENDIF … … 323 323 ALLOCATE( sf_t_rnf(1)%fnow(jpi,jpj,1) ) 324 324 IF( sn_t_rnf%ln_tint ) ALLOCATE( sf_t_rnf(1)%fdta(jpi,jpj,1,2) ) 325 CALL fld_fill (sf_t_rnf, (/ sn_t_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf' ) 325 CALL fld_fill (sf_t_rnf, (/ sn_t_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf' ) 326 326 ENDIF 327 327 ! … … 335 335 ALLOCATE( sf_s_rnf(1)%fnow(jpi,jpj,1) ) 336 336 IF( sn_s_rnf%ln_tint ) ALLOCATE( sf_s_rnf(1)%fdta(jpi,jpj,1,2) ) 337 CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf' ) 338 ENDIF 339 ! 340 IF( ln_rnf_depth ) THEN ! depth of runoffs set from a file 337 CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf' ) 338 ENDIF 339 ! 340 IF( ln_rnf_depth ) THEN ! depth of runoffs set from a file 341 341 IF(lwp) WRITE(numout,*) 342 342 IF(lwp) WRITE(numout,*) ' runoffs depth read in a file' 343 rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname ) 344 CALL iom_open ( rn_dep_file, inum ) ! open file 345 CALL iom_get ( inum, jpdom_data, sn_dep_rnf%clvar, h_rnf ) ! read the river mouth array 346 CALL iom_close( inum ) ! close file 343 rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname ) 344 CALL iom_open ( rn_dep_file, inum ) ! open file 345 CALL iom_get ( inum, jpdom_data, sn_dep_rnf%clvar, h_rnf ) ! read the river mouth array 346 CALL iom_close( inum ) ! close file 347 347 ! 348 348 nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied 349 DO jj = 1, jpj 350 DO ji = 1, jpi 351 IF( h_rnf(ji,jj) > 0._wp ) THEN 352 jk = 2 353 DO WHILE ( jk /= mbkt(ji,jj) .AND. fsdept(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1 ; END DO 354 nk_rnf(ji,jj) = jk 355 ELSEIF( h_rnf(ji,jj) == -1 ) THEN ; nk_rnf(ji,jj) = 1 349 DO jj = 1, jpj 350 DO ji = 1, jpi 351 IF( h_rnf(ji,jj) > 0._wp ) THEN 352 jk = 2 353 DO WHILE ( jk /= mbkt(ji,jj) .AND. fsdept(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1 ; END DO 354 nk_rnf(ji,jj) = jk 355 ELSEIF( h_rnf(ji,jj) == -1 ) THEN ; nk_rnf(ji,jj) = 1 356 356 ELSEIF( h_rnf(ji,jj) == -999 ) THEN ; nk_rnf(ji,jj) = mbkt(ji,jj) 357 ELSEIF( h_rnf(ji,jj) /= 0 ) THEN 358 CALL ctl_stop( 'runoff depth not positive, and not -999 or -1, rnf value in file fort.999' ) 359 WRITE(999,*) 'ji, jj, rnf(ji,jj) :', ji, jj, rnf(ji,jj) 360 ENDIF 361 END DO 362 END DO 363 DO jj = 1, jpj ! set the associated depth 364 DO ji = 1, jpi 357 ELSEIF( h_rnf(ji,jj) /= 0 ) THEN 358 CALL ctl_stop( 'runoff depth not positive, and not -999 or -1, rnf value in file fort.999' ) 359 WRITE(999,*) 'ji, jj, rnf(ji,jj) :', ji, jj, rnf(ji,jj) 360 ENDIF 361 END DO 362 END DO 363 DO jj = 1, jpj ! set the associated depth 364 DO ji = 1, jpi 365 365 h_rnf(ji,jj) = 0._wp 366 DO jk = 1, nk_rnf(ji,jj) 367 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) 366 DO jk = 1, nk_rnf(ji,jj) 367 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) 368 368 END DO 369 369 END DO 370 370 END DO 371 ELSE ! runoffs applied at the surface 372 nk_rnf(:,:) = 1 371 ELSE ! runoffs applied at the surface 372 nk_rnf(:,:) = 1 373 373 h_rnf (:,:) = fse3t(:,:,1) 374 ENDIF 375 ! 374 ENDIF 375 ! 376 376 ENDIF 377 377 ! … … 389 389 ! 390 390 IF ( ln_rnf_depth ) CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already', & 391 & 'be spread through depth by ln_rnf_depth' ) 391 & 'be spread through depth by ln_rnf_depth' ) 392 392 ! 393 393 nkrnf = 0 ! Number of level over which Kz increase … … 410 410 IF(lwp) WRITE(numout,*) 411 411 IF(lwp) WRITE(numout,*) ' No specific treatment at river mouths' 412 rnfmsk (:,:) = 0._wp 412 rnfmsk (:,:) = 0._wp 413 413 rnfmsk_z(:) = 0._wp 414 414 nkrnf = 0 … … 421 421 !!---------------------------------------------------------------------- 422 422 !! *** ROUTINE rnf_mouth *** 423 !! 423 !! 424 424 !! ** Purpose : define the river mouths mask 425 425 !! 426 426 !! ** Method : read the river mouth mask (=0/1) in the river runoff 427 !! climatological file. Defined a given vertical structure. 428 !! CAUTION, the vertical structure is hard coded on the 427 !! climatological file. Defined a given vertical structure. 428 !! CAUTION, the vertical structure is hard coded on the 429 429 !! first 5 levels. 430 430 !! This fields can be used to: 431 !! - set an upstream advection scheme 431 !! - set an upstream advection scheme 432 432 !! (ln_rnf_mouth=T and ln_traadv_cen2=T) 433 !! - increase vertical on the top nn_krnf vertical levels 433 !! - increase vertical on the top nn_krnf vertical levels 434 434 !! at river runoff input grid point (nn_krnf>=2, see step.F90) 435 435 !! - set to zero SSS restoring flux at river mouth grid points … … 442 442 CHARACTER(len=140) :: cl_rnfile ! runoff file name 443 443 !!---------------------------------------------------------------------- 444 ! 444 ! 445 445 IF(lwp) WRITE(numout,*) 446 446 IF(lwp) WRITE(numout,*) 'rnf_mouth : river mouth mask' … … 451 451 IF( sn_cnf%cltype == 'monthly' ) WRITE(cl_rnfile, '(a,"m",i2)' ) TRIM( cl_rnfile ), nmonth ! add month 452 452 ENDIF 453 453 454 454 ! horizontal mask (read in NetCDF file) 455 455 CALL iom_open ( cl_rnfile, inum ) ! open file 456 456 CALL iom_get ( inum, jpdom_data, sn_cnf%clvar, rnfmsk ) ! read the river mouth array 457 457 CALL iom_close( inum ) ! close file 458 458 459 459 IF( nn_closea == 1 ) CALL clo_rnf( rnfmsk ) ! closed sea inflow set as ruver mouth 460 460 461 rnfmsk_z(:) = 0._wp ! vertical structure 461 rnfmsk_z(:) = 0._wp ! vertical structure 462 462 rnfmsk_z(1) = 1.0 463 463 rnfmsk_z(2) = 1.0 ! ********** … … 465 465 rnfmsk_z(4) = 0.25 ! ********** 466 466 rnfmsk_z(5) = 0.125 467 ! 467 ! 468 468 END SUBROUTINE rnf_mouth 469 469 470 470 !!====================================================================== 471 471 END MODULE sbcrnf -
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
r3294 r3598 8 8 !! NEMO 1.0 ! 2002-08 (G. Madec) free form + modules 9 9 !! - ! 2004-01 (A. de Miranda, G. Madec, J.M. Molines ) add advective bbl 10 !! 3.3 ! 2009-11 (G. Madec) merge trabbl and trabbl_adv + style + optimization 11 !! - ! 2010-04 (G. Madec) Campin & Goosse advective bbl 10 !! 3.3 ! 2009-11 (G. Madec) merge trabbl and trabbl_adv + style + optimization 11 !! - ! 2010-04 (G. Madec) Campin & Goosse advective bbl 12 12 !! - ! 2010-06 (C. Ethe, G. Madec) merge TRA-TRC 13 13 !! - ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level … … 30 30 USE trdmod_oce ! trends: ocean variables 31 31 USE trdtra ! trends: active tracers 32 USE iom ! IOM server 32 USE iom ! IOM server 33 33 USE in_out_manager ! I/O manager 34 34 USE lbclnk ! ocean lateral boundary conditions … … 49 49 LOGICAL, PUBLIC, PARAMETER :: lk_trabbl = .TRUE. !: bottom boundary layer flag 50 50 51 ! !!* Namelist nambbl * 51 ! !!* Namelist nambbl * 52 52 INTEGER , PUBLIC :: nn_bbl_ldf = 0 !: =1 : diffusive bbl or not (=0) 53 53 INTEGER , PUBLIC :: nn_bbl_adv = 0 !: =1/2 : advective bbl or not (=0) … … 58 58 59 59 LOGICAL , PUBLIC :: l_bbl !: flag to compute bbl diffu. flux coef and transport 60 60 61 61 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: utr_bbl , vtr_bbl ! u- (v-) transport in the bottom boundary layer 62 62 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: ahu_bbl , ahv_bbl ! masked diffusive bbl coeff. at u & v-pts 63 63 64 INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku_d , mbkv_d ! vertical index of the "lower" bottom ocean U/V-level65 INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: mgrhu , mgrhv ! = +/-1, sign of grad(H) in u-(v-)direction66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ahu_bbl_0, ahv_bbl_0 ! diffusive bbl flux coefficients at u and v-points67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3u_bbl_0, e3v_bbl_0 ! thichness of the bbl (e3) at u and v-points68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_e1e2t ! inverse of the cell surface at t-point [1/m2]64 INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: mbku_d , mbkv_d ! vertical index of the "lower" bottom ocean U/V-level (PUBLIC for TAM) 65 INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: mgrhu , mgrhv ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM) 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ahu_bbl_0, ahv_bbl_0 ! diffusive bbl flux coefficients at u and v-points 67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: e3u_bbl_0, e3v_bbl_0 ! thichness of the bbl (e3) at u and v-points (PUBLIC for TAM) 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: r1_e1e2t ! inverse of the cell surface at t-point [1/m2] (PUBLIC for TAM) 69 69 70 70 !! * Substitutions … … 95 95 !!---------------------------------------------------------------------- 96 96 !! *** ROUTINE bbl *** 97 !! 98 !! ** Purpose : Compute the before tracer (t & s) trend associated 97 !! 98 !! ** Purpose : Compute the before tracer (t & s) trend associated 99 99 !! with the bottom boundary layer and add it to the general 100 100 !! trend of tracer equations. … … 103 103 !! diffusive and/or advective contribution to the tracer trend 104 104 !! is added to the general tracer trend 105 !!---------------------------------------------------------------------- 106 INTEGER, INTENT( in ) :: kt ! ocean time-step 105 !!---------------------------------------------------------------------- 106 INTEGER, INTENT( in ) :: kt ! ocean time-step 107 107 !! 108 108 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds … … 112 112 ! 113 113 IF( l_trdtra ) THEN !* Save ta and sa trends 114 CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 115 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 114 CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 115 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 116 116 ztrds(:,:,:) = tsa(:,:,:,jp_sal) 117 117 ENDIF 118 118 119 119 IF( l_bbl ) CALL bbl( kt, nit000, 'TRA' ) !* bbl coef. and transport (only if not already done in trcbbl) 120 120 121 121 IF( nn_bbl_ldf == 1 ) THEN !* Diffusive bbl 122 122 ! … … 125 125 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf - Ta: ', mask1=tmask, & 126 126 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 127 ! lateral boundary conditions ; just need for outputs 127 ! lateral boundary conditions ; just need for outputs 128 128 CALL lbc_lnk( ahu_bbl, 'U', 1. ) ; CALL lbc_lnk( ahv_bbl, 'V', 1. ) 129 CALL iom_put( "ahu_bbl", ahu_bbl ) ! bbl diffusive flux i-coef 129 CALL iom_put( "ahu_bbl", ahu_bbl ) ! bbl diffusive flux i-coef 130 130 CALL iom_put( "ahv_bbl", ahv_bbl ) ! bbl diffusive flux j-coef 131 131 ! … … 138 138 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv - Ta: ', mask1=tmask, & 139 139 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 140 ! lateral boundary conditions ; just need for outputs 140 ! lateral boundary conditions ; just need for outputs 141 141 CALL lbc_lnk( utr_bbl, 'U', 1. ) ; CALL lbc_lnk( vtr_bbl, 'V', 1. ) 142 CALL iom_put( "uoce_bbl", utr_bbl ) ! bbl i-transport 142 CALL iom_put( "uoce_bbl", utr_bbl ) ! bbl i-transport 143 143 CALL iom_put( "voce_bbl", vtr_bbl ) ! bbl j-transport 144 144 ! … … 150 150 CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbl, ztrdt ) 151 151 CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_bbl, ztrds ) 152 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 152 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 153 153 ENDIF 154 154 ! … … 161 161 !!---------------------------------------------------------------------- 162 162 !! *** ROUTINE tra_bbl_dif *** 163 !! 163 !! 164 164 !! ** Purpose : Computes the bottom boundary horizontal and vertical 165 !! advection terms. 166 !! 167 !! ** Method : 165 !! advection terms. 166 !! 167 !! ** Method : 168 168 !! * diffusive bbl (nn_bbl_ldf=1) : 169 169 !! When the product grad( rho) * grad(h) < 0 (where grad is an … … 179 179 !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. 180 180 !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 181 !!---------------------------------------------------------------------- 181 !!---------------------------------------------------------------------- 182 182 ! 183 183 INTEGER , INTENT(in ) :: kjpt ! number of tracers 184 184 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 185 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 185 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 186 186 ! 187 187 INTEGER :: ji, jj, jn ! dummy loop indices … … 202 202 #else 203 203 DO jj = 1, jpj 204 DO ji = 1, jpi 204 DO ji = 1, jpi 205 205 #endif 206 206 ik = mbkt(ji,jj) ! bottom T-level index … … 233 233 ! 234 234 END SUBROUTINE tra_bbl_dif 235 235 236 236 237 237 SUBROUTINE tra_bbl_adv( ptb, pta, kjpt ) … … 239 239 !! *** ROUTINE trc_bbl *** 240 240 !! 241 !! ** Purpose : Compute the before passive tracer trend associated 241 !! ** Purpose : Compute the before passive tracer trend associated 242 242 !! with the bottom boundary layer and add it to the general trend 243 243 !! of tracer equations. 244 244 !! ** Method : advective bbl (nn_bbl_adv = 1 or 2) : 245 245 !! nn_bbl_adv = 1 use of the ocean near bottom velocity as bbl velocity 246 !! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation i.e. 247 !! transport proportional to the along-slope density gradient 246 !! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation i.e. 247 !! transport proportional to the along-slope density gradient 248 248 !! 249 249 !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. 250 250 !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 251 !!---------------------------------------------------------------------- 251 !!---------------------------------------------------------------------- 252 252 INTEGER , INTENT(in ) :: kjpt ! number of tracers 253 253 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 254 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 254 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 255 255 ! 256 256 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 264 264 ! ! =========== 265 265 DO jn = 1, kjpt ! tracer loop 266 ! ! =========== 266 ! ! =========== 267 267 # if defined key_vectopt_loop 268 268 DO jj = 1, 1 … … 282 282 ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr 283 283 pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra 284 ! 284 ! 285 285 DO jk = ikus, ikud-1 ! down-slope upper to down T-point (deep column) 286 286 zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,jk) … … 288 288 pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra 289 289 END DO 290 ! 290 ! 291 291 zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,ikud) 292 292 ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr … … 299 299 ikvd = mbkv_d(ji,jj) ; ikvs = mbkv(ji,jj) 300 300 zv_bbl = ABS( vtr_bbl(ji,jj) ) 301 ! 301 ! 302 302 ! up -slope T-point (shelf bottom point) 303 303 zbtr = r1_e1e2t(ji,ijs) / fse3t(ji,ijs,ikvs) 304 304 ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr 305 305 pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra 306 ! 306 ! 307 307 DO jk = ikvs, ikvd-1 ! down-slope upper to down T-point (deep column) 308 308 zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,jk) … … 330 330 !!---------------------------------------------------------------------- 331 331 !! *** ROUTINE bbl *** 332 !! 332 !! 333 333 !! ** Purpose : Computes the bottom boundary horizontal and vertical 334 !! advection terms. 335 !! 336 !! ** Method : 334 !! advection terms. 335 !! 336 !! ** Method : 337 337 !! * diffusive bbl (nn_bbl_ldf=1) : 338 338 !! When the product grad( rho) * grad(h) < 0 (where grad is an … … 353 353 !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. 354 354 !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 355 !!---------------------------------------------------------------------- 355 !!---------------------------------------------------------------------- 356 356 ! 357 357 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 399 399 - 0.121555e-07 ) * pfh 400 400 !!---------------------------------------------------------------------- 401 401 402 402 ! 403 403 IF( nn_timing == 1 ) CALL timing_start( 'bbl') 404 404 ! 405 CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 406 ! 407 405 CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 406 ! 407 408 408 IF( kt == kit000 ) THEN 409 409 IF(lwp) WRITE(numout,*) … … 411 411 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 412 412 ENDIF 413 413 414 414 ! !* bottom temperature, salinity, velocity and depth 415 415 #if defined key_vectopt_loop … … 426 426 ! 427 427 zub(ji,jj) = un(ji,jj,mbku(ji,jj)) ! bottom velocity 428 zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) 428 zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) 429 429 END DO 430 430 END DO 431 431 432 432 ! !-------------------! 433 433 IF( nn_bbl_ldf == 1 ) THEN ! diffusive bbl ! 434 434 ! !-------------------! 435 435 DO jj = 1, jpjm1 ! (criteria for non zero flux: grad(rho).grad(h) < 0 ) 436 DO ji = 1, jpim1 437 ! ! i-direction 436 DO ji = 1, jpim1 437 ! ! i-direction 438 438 zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) ) ! T, S anomalie, and depth 439 439 zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 … … 442 442 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) ) & 443 443 & - ( zsb(ji+1,jj) - zsb(ji,jj) ) ) * umask(ji,jj,1) 444 ! 444 ! 445 445 zsign = SIGN( 0.5, - zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of ( i-gradient * i-slope ) 446 446 ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj) ! masked diffusive flux coeff. 447 447 ! 448 ! ! j-direction 448 ! ! j-direction 449 449 zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) ) ! T, S anomalie, and depth 450 450 zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0 … … 453 453 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) ) & 454 454 & - ( zsb(ji,jj+1) - zsb(ji,jj) ) ) * vmask(ji,jj,1) 455 ! 455 ! 456 456 zsign = SIGN( 0.5, -zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of ( j-gradient * j-slope ) 457 457 ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj) … … 475 475 zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 476 476 ! ! masked bbl i-gradient of density 477 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) ) & 477 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) ) & 478 478 & - ( zsb(ji+1,jj) - zsb(ji,jj) ) ) * umask(ji,jj,1) 479 ! 479 ! 480 480 zsign = SIGN( 0.5, - zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of i-gradient * i-slope 481 481 zsigna= SIGN( 0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) ) ) ! sign of u * i-slope … … 489 489 zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 490 490 ! ! masked bbl j-gradient of density 491 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) ) & 491 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) ) & 492 492 & - ( zsb(ji,jj+1) - zsb(ji,jj) ) ) * vmask(ji,jj,1) 493 493 zsign = SIGN( 0.5, - zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of j-gradient * j-slope … … 513 513 zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 514 514 zgdrho = fsbeta( zt, zs, zh ) & 515 & * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) ) & 515 & * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) ) & 516 516 & - ( zsb(iid,jj) - zsb(iis,jj) ) ) * umask(ji,jj,1) 517 517 zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep … … 530 530 zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) ) 531 531 zgdrho = fsbeta( zt, zs, zh ) & 532 & * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) ) & 532 & * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) ) & 533 533 & - ( zsb(ji,ijd) - zsb(ji,ijs) ) ) * vmask(ji,jj,1) 534 534 zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep … … 542 542 ENDIF 543 543 ! 544 CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 544 CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 545 545 ! 546 546 IF( nn_timing == 1 ) CALL timing_stop( 'bbl') … … 567 567 IF( nn_timing == 1 ) CALL timing_start( 'tra_bbl_init') 568 568 ! 569 CALL wrk_alloc( jpi, jpj, zmbk ) 569 CALL wrk_alloc( jpi, jpj, zmbk ) 570 570 ! 571 571 … … 588 588 ! ! allocate trabbl arrays 589 589 IF( tra_bbl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' ) 590 590 591 591 IF( nn_bbl_adv == 1 ) WRITE(numout,*) ' * Advective BBL using upper velocity' 592 592 IF( nn_bbl_adv == 2 ) WRITE(numout,*) ' * Advective BBL using velocity = F( delta rho)' … … 597 597 ! !* inverse of surface of T-cells 598 598 r1_e1e2t(:,:) = 1._wp / ( e1t(:,:) * e2t(:,:) ) 599 599 600 600 ! !* vertical index of "deep" bottom u- and v-points 601 601 DO jj = 1, jpjm1 ! (the "shelf" bottom k-indices are mbku and mbkv) … … 605 605 END DO 606 606 END DO 607 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 607 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 608 608 zmbk(:,:) = REAL( mbku_d(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 609 609 zmbk(:,:) = REAL( mbkv_d(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) … … 611 611 !* sign of grad(H) at u- and v-points 612 612 mgrhu(jpi,:) = 0. ; mgrhu(:,jpj) = 0. ; mgrhv(jpi,:) = 0. ; mgrhv(:,jpj) = 0. 613 DO jj = 1, jpjm1 613 DO jj = 1, jpjm1 614 614 DO ji = 1, jpim1 615 615 mgrhu(ji,jj) = INT( SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) ) ) … … 618 618 END DO 619 619 620 DO jj = 1, jpjm1 !* bbl thickness at u- (v-) point 620 DO jj = 1, jpjm1 !* bbl thickness at u- (v-) point 621 621 DO ji = 1, jpim1 ! minimum of top & bottom e3u_0 (e3v_0) 622 e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj )), fse3u_0(ji,jj,mbkt(ji,jj)) ) 623 e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) ) 624 END DO 622 e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj )), fse3u_0(ji,jj,mbkt(ji,jj)) ) 623 e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) ) 624 END DO 625 625 END DO 626 626 CALL lbc_lnk( e3u_bbl_0, 'U', 1. ) ; CALL lbc_lnk( e3v_bbl_0, 'V', 1. ) ! lateral boundary conditions 627 627 628 ! !* masked diffusive flux coefficients 628 ! !* masked diffusive flux coefficients 629 629 ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:) * umask(:,:,1) 630 630 ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:) * vmask(:,:,1) … … 636 636 CASE ( 2 ) ! ORCA_R2 637 637 ij0 = 102 ; ij1 = 102 ! Gibraltar enhancement of BBL 638 ii0 = 139 ; ii1 = 140 638 ii0 = 139 ; ii1 = 140 639 639 ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) 640 640 ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) … … 647 647 CASE ( 4 ) ! ORCA_R4 648 648 ij0 = 52 ; ij1 = 52 ! Gibraltar enhancement of BBL 649 ii0 = 70 ; ii1 = 71 649 ii0 = 70 ; ii1 = 71 650 650 ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) 651 651 ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) … … 654 654 ENDIF 655 655 ! 656 CALL wrk_dealloc( jpi, jpj, zmbk ) 656 CALL wrk_dealloc( jpi, jpj, zmbk ) 657 657 ! 658 658 IF( nn_timing == 1 ) CALL timing_stop( 'tra_bbl_init') -
trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r3294 r3598 16 16 !!---------------------------------------------------------------------- 17 17 USE oce ! ocean dynamics and tracers variables 18 USE dom_oce ! ocean space and time domain variables 18 USE dom_oce ! ocean space and time domain variables 19 19 USE zdf_oce ! ocean vertical physics variables 20 20 USE in_out_manager ! I/O manager … … 31 31 32 32 ! !!* Namelist nambfr: bottom friction namelist * 33 INTEGER :: nn_bfr = 0 ! = 0/1/2/3 type of bottom friction34 REAL(wp) :: rn_bfri1 = 4.0e-4_wp ! bottom drag coefficient (linear case)35 REAL(wp) :: rn_bfri2 = 1.0e-3_wp ! bottom drag coefficient (non linear case)36 REAL(wp) :: rn_bfeb2 = 2.5e-3_wp ! background bottom turbulent kinetic energy [m2/s2]37 REAL(wp) :: rn_bfrien = 30._wp ! local factor to enhance coefficient bfri38 LOGICAL :: ln_bfr2d = .false. ! logical switch for 2D enhancement39 LOGICAL , PUBLIC :: ln_bfrimp = .false. ! logical switch for implicit bottom friction40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bfrcoef2d ! 2D bottom drag coefficient33 INTEGER , PUBLIC :: nn_bfr = 0 ! = 0/1/2/3 type of bottom friction (PUBLIC for TAM) 34 REAL(wp), PUBLIC :: rn_bfri1 = 4.0e-4_wp ! bottom drag coefficient (linear case) (PUBLIC for TAM) 35 REAL(wp), PUBLIC :: rn_bfri2 = 1.0e-3_wp ! bottom drag coefficient (non linear case) (PUBLIC for TAM) 36 REAL(wp), PUBLIC :: rn_bfeb2 = 2.5e-3_wp ! background bottom turbulent kinetic energy [m2/s2] (PUBLIC for TAM) 37 REAL(wp), PUBLIC :: rn_bfrien = 30._wp ! local factor to enhance coefficient bfri (PUBLIC for TAM) 38 LOGICAL , PUBLIC :: ln_bfr2d = .false. ! logical switch for 2D enhancement (PUBLIC for TAM) 39 LOGICAL , PUBLIC :: ln_bfrimp = .false. ! logical switch for implicit bottom friction 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: bfrcoef2d ! 2D bottom drag coefficient (PUBLIC for TAM) 41 41 42 42 !! * Substitutions … … 64 64 !!---------------------------------------------------------------------- 65 65 !! *** ROUTINE zdf_bfr *** 66 !! 66 !! 67 67 !! ** Purpose : compute the bottom friction coefficient. 68 68 !! 69 !! ** Method : Calculate and store part of the momentum trend due 70 !! to bottom friction following the chosen friction type 69 !! ** Method : Calculate and store part of the momentum trend due 70 !! to bottom friction following the chosen friction type 71 71 !! (free-slip, linear, or quadratic). The component 72 72 !! calculated here is multiplied by the bottom velocity in … … 102 102 DO ji = 2, jpim1 103 103 # endif 104 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 104 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 105 105 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 106 106 ! … … 113 113 zecv = SQRT( vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 114 114 ! 115 bfrua(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj ) ) * zecu 115 bfrua(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj ) ) * zecu 116 116 bfrva(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji ,jj+1) ) * zecv 117 117 END DO … … 133 133 !!---------------------------------------------------------------------- 134 134 !! *** ROUTINE zdf_bfr_init *** 135 !! 135 !! 136 136 !! ** Purpose : Initialization of the bottom friction 137 137 !! … … 193 193 bfrcoef2d(:,:) = rn_bfri1 ! initialize bfrcoef2d to the namelist variable 194 194 ! 195 IF(ln_bfr2d) THEN 195 IF(ln_bfr2d) THEN 196 196 ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 197 197 CALL iom_open('bfr_coef.nc',inum) … … 213 213 bfrcoef2d(:,:) = rn_bfri2 ! initialize bfrcoef2d to the namelist variable 214 214 ! 215 IF(ln_bfr2d) THEN 215 IF(ln_bfr2d) THEN 216 216 ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 217 217 CALL iom_open('bfr_coef.nc',inum)
Note: See TracChangeset
for help on using the changeset viewer.
