Changeset 3359
 Timestamp:
 Oct 16, 2018 7:36:26 PM (6 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/turbulence_closure_mod.f90
r3300 r3359 25 25 !  26 26 ! $Id$ 27 ! Restructured loops and ifs in production_e to ease vectorization on GPUs 28 ! 29 ! 3300 20181002 14:16:54Z gronemeier 27 30 !  removed global array wall_flags_0_global, hence reduced accuracy of l_wall 28 31 ! calculation … … 2570 2573 INTEGER(iwp) :: surf_e !< end index of surface elements at given ij position 2571 2574 INTEGER(iwp) :: surf_s !< start index of surface elements at given ij position 2575 INTEGER(iwp) :: flag_nr !< number of masking flag 2572 2576 2573 2577 REAL(wp) :: def !< … … 2584 2588 REAL(wp) :: wsvs !< momentum flux w"v" 2585 2589 2586 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudx !< Gradient of ucomponent in xdirection 2587 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudy !< Gradient of ucomponent in ydirection 2588 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudz !< Gradient of ucomponent in zdirection 2589 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdx !< Gradient of vcomponent in xdirection 2590 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdy !< Gradient of vcomponent in ydirection 2591 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdz !< Gradient of vcomponent in zdirection 2592 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdx !< Gradient of wcomponent in xdirection 2593 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdy !< Gradient of wcomponent in ydirection 2594 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdz !< Gradient of wcomponent in zdirection 2595 2590 REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of ucomponent in xdirection 2591 REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of ucomponent in ydirection 2592 REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of ucomponent in zdirection 2593 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of vcomponent in xdirection 2594 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of vcomponent in ydirection 2595 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of vcomponent in zdirection 2596 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of wcomponent in xdirection 2597 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of wcomponent in ydirection 2598 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of wcomponent in zdirection 2599 2600 2601 2602 ! 2603 ! Calculate TKE production by shear. Calculate gradients at all grid 2604 ! points first, gradients at surfacebounded grid points will be 2605 ! overwritten further below. 2596 2606 DO i = nxl, nxr 2597 2598 IF ( constant_flux_layer ) THEN 2599 2600 ! 2601 ! Calculate TKE production by shear. Calculate gradients at all grid 2602 ! points first, gradients at surfacebounded grid points will be 2603 ! overwritten further below. 2604 DO j = nys, nyn 2605 DO k = nzb+1, nzt 2606 2607 dudx(k,j) = ( u(k,j,i+1)  u(k,j,i) ) * ddx 2608 dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1)  & 2609 u(k,j1,i)  u(k,j1,i+1) ) * ddy 2610 dudz(k,j) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)  & 2611 u(k1,j,i)  u(k1,j,i+1) ) * & 2612 dd2zu(k) 2613 2614 dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1)  & 2615 v(k,j,i1)  v(k,j+1,i1) ) * ddx 2616 dvdy(k,j) = ( v(k,j+1,i)  v(k,j,i) ) * ddy 2617 dvdz(k,j) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)  & 2618 v(k1,j,i)  v(k1,j+1,i) ) * & 2619 dd2zu(k) 2620 2621 dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k1,j,i+1)  & 2622 w(k,j,i1)  w(k1,j,i1) ) * ddx 2623 dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k1,j+1,i)  & 2624 w(k,j1,i)  w(k1,j1,i) ) * ddy 2625 dwdz(k,j) = ( w(k,j,i)  w(k1,j,i) ) * ddzw(k) 2626 2627 ENDDO 2607 DO j = nys, nyn 2608 DO k = nzb+1, nzt 2609 2610 dudx(k) = ( u(k,j,i+1)  u(k,j,i) ) * ddx 2611 dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1)  & 2612 u(k,j1,i)  u(k,j1,i+1) ) * ddy 2613 dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)  & 2614 u(k1,j,i)  u(k1,j,i+1) ) * dd2zu(k) 2615 2616 dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1)  & 2617 v(k,j,i1)  v(k,j+1,i1) ) * ddx 2618 dvdy(k) = ( v(k,j+1,i)  v(k,j,i) ) * ddy 2619 dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)  & 2620 v(k1,j,i)  v(k1,j+1,i) ) * dd2zu(k) 2621 2622 dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k1,j,i+1)  & 2623 w(k,j,i1)  w(k1,j,i1) ) * ddx 2624 dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k1,j+1,i)  & 2625 w(k,j1,i)  w(k1,j1,i) ) * ddy 2626 dwdz(k) = ( w(k,j,i)  w(k1,j,i) ) * ddzw(k) 2627 2628 2628 ENDDO 2629 2629 2630 ! 2631 ! Position beneath wall 2632 ! (2)  Will allways be executed. 2633 ! 'bottom and wall: use u_0,v_0 and wall functions' 2634 DO j = nys, nyn 2630 2631 flag_nr = 29 2632 2633 2634 IF ( constant_flux_layer ) THEN 2635 ! 2636 2637 flag_nr = 0 2638 2639 ! Position beneath wall 2640 ! (2)  Will allways be executed. 2641 ! 'bottom and wall: use u_0,v_0 and wall functions' 2635 2642 ! 2636 2643 ! Compute gradients at north and southfacing surfaces. 2637 ! First, for default surfaces, then for urban surfaces. 2644 ! First, for default surfaces, then for urban surfaces. 2638 2645 ! Note, so far no natural vertical surfaces implemented 2639 2646 DO l = 0, 1 … … 2644 2651 usvs = surf_def_v(l)%mom_flux_tke(0,m) 2645 2652 wsvs = surf_def_v(l)%mom_flux_tke(1,m) 2646 2653 2647 2654 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 2648 2655 * 0.5_wp * dy … … 2650 2657 ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 2651 2658 sign_dir = MERGE( 1.0_wp, 1.0_wp, & 2652 BTEST( wall_flags_0(k,j1,i), 0 ) )2653 dudy(k ,j) = sign_dir * usvs / ( km_neutral + 1E10_wp )2654 dwdy(k ,j) = sign_dir * wsvs / ( km_neutral + 1E10_wp )2659 BTEST( wall_flags_0(k,j1,i), flag_nr ) ) 2660 dudy(k) = sign_dir * usvs / ( km_neutral + 1E10_wp ) 2661 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E10_wp ) 2655 2662 ENDDO 2656 2663 ! … … 2662 2669 usvs = surf_lsm_v(l)%mom_flux_tke(0,m) 2663 2670 wsvs = surf_lsm_v(l)%mom_flux_tke(1,m) 2664 2671 2665 2672 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 2666 2673 * 0.5_wp * dy … … 2668 2675 ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 2669 2676 sign_dir = MERGE( 1.0_wp, 1.0_wp, & 2670 BTEST( wall_flags_0(k,j1,i), 0 ) )2671 dudy(k ,j) = sign_dir * usvs / ( km_neutral + 1E10_wp )2672 dwdy(k ,j) = sign_dir * wsvs / ( km_neutral + 1E10_wp )2673 ENDDO 2677 BTEST( wall_flags_0(k,j1,i), flag_nr ) ) 2678 dudy(k) = sign_dir * usvs / ( km_neutral + 1E10_wp ) 2679 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E10_wp ) 2680 ENDDO 2674 2681 ! 2675 2682 ! Urban surfaces … … 2680 2687 usvs = surf_usm_v(l)%mom_flux_tke(0,m) 2681 2688 wsvs = surf_usm_v(l)%mom_flux_tke(1,m) 2682 2689 2683 2690 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 2684 2691 * 0.5_wp * dy … … 2686 2693 ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 2687 2694 sign_dir = MERGE( 1.0_wp, 1.0_wp, & 2688 BTEST( wall_flags_0(k,j1,i), 0 ) )2689 dudy(k ,j) = sign_dir * usvs / ( km_neutral + 1E10_wp )2690 dwdy(k ,j) = sign_dir * wsvs / ( km_neutral + 1E10_wp )2691 ENDDO 2695 BTEST( wall_flags_0(k,j1,i), flag_nr ) ) 2696 dudy(k) = sign_dir * usvs / ( km_neutral + 1E10_wp ) 2697 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E10_wp ) 2698 ENDDO 2692 2699 ENDDO 2693 2700 ! … … 2706 2713 ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 2707 2714 sign_dir = MERGE( 1.0_wp, 1.0_wp, & 2708 BTEST( wall_flags_0(k,j,i1), 0 ) )2709 dvdx(k ,j) = sign_dir * vsus / ( km_neutral + 1E10_wp )2710 dwdx(k ,j) = sign_dir * wsus / ( km_neutral + 1E10_wp )2711 ENDDO 2712 ! 2713 ! Natural surfaces 2715 BTEST( wall_flags_0(k,j,i1), flag_nr ) ) 2716 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E10_wp ) 2717 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E10_wp ) 2718 ENDDO 2719 ! 2720 ! Natural surfaces 2714 2721 surf_s = surf_lsm_v(l)%start_index(j,i) 2715 2722 surf_e = surf_lsm_v(l)%end_index(j,i) … … 2724 2731 ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 2725 2732 sign_dir = MERGE( 1.0_wp, 1.0_wp, & 2726 BTEST( wall_flags_0(k,j,i1), 0 ) )2727 dvdx(k ,j) = sign_dir * vsus / ( km_neutral + 1E10_wp )2728 dwdx(k ,j) = sign_dir * wsus / ( km_neutral + 1E10_wp )2729 ENDDO 2730 ! 2731 ! Urban surfaces 2733 BTEST( wall_flags_0(k,j,i1), flag_nr ) ) 2734 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E10_wp ) 2735 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E10_wp ) 2736 ENDDO 2737 ! 2738 ! Urban surfaces 2732 2739 surf_s = surf_usm_v(l)%start_index(j,i) 2733 2740 surf_e = surf_usm_v(l)%end_index(j,i) … … 2742 2749 ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 2743 2750 sign_dir = MERGE( 1.0_wp, 1.0_wp, & 2744 BTEST( wall_flags_0(k,j,i1), 0 ) )2745 dvdx(k ,j) = sign_dir * vsus / ( km_neutral + 1E10_wp )2746 dwdx(k ,j) = sign_dir * wsus / ( km_neutral + 1E10_wp )2747 ENDDO 2751 BTEST( wall_flags_0(k,j,i1), flag_nr ) ) 2752 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E10_wp ) 2753 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E10_wp ) 2754 ENDDO 2748 2755 ENDDO 2749 2756 ! … … 2754 2761 k = surf_def_h(0)%k(m) 2755 2762 ! 2756 ! Please note, actually, an interpolation of u_0 and v_0 2757 ! onto the grid center would be required. However, this 2763 ! Please note, actually, an interpolation of u_0 and v_0 2764 ! onto the grid center would be required. However, this 2758 2765 ! would require several data transfers between 2Dgrid and 2759 ! wall type. The effect of this missing interpolation is 2766 ! wall type. The effect of this missing interpolation is 2760 2767 ! negligible. (See also production_e_init). 2761 dudz(k ,j) = ( u(k+1,j,i)  surf_def_h(0)%u_0(m) ) * dd2zu(k)2762 dvdz(k ,j) = ( v(k+1,j,i)  surf_def_h(0)%v_0(m) ) * dd2zu(k)2763 2768 dudz(k) = ( u(k+1,j,i)  surf_def_h(0)%u_0(m) ) * dd2zu(k) 2769 dvdz(k) = ( v(k+1,j,i)  surf_def_h(0)%v_0(m) ) * dd2zu(k) 2770 2764 2771 ENDDO 2765 2772 ! … … 2770 2777 k = surf_lsm_h%k(m) 2771 2778 2772 dudz(k ,j) = ( u(k+1,j,i)  surf_lsm_h%u_0(m) ) * dd2zu(k)2773 dvdz(k ,j) = ( v(k+1,j,i)  surf_lsm_h%v_0(m) ) * dd2zu(k)2774 2779 dudz(k) = ( u(k+1,j,i)  surf_lsm_h%u_0(m) ) * dd2zu(k) 2780 dvdz(k) = ( v(k+1,j,i)  surf_lsm_h%v_0(m) ) * dd2zu(k) 2781 2775 2782 ENDDO 2776 2783 ! … … 2781 2788 k = surf_usm_h%k(m) 2782 2789 2783 dudz(k ,j) = ( u(k+1,j,i)  surf_usm_h%u_0(m) ) * dd2zu(k)2784 dvdz(k ,j) = ( v(k+1,j,i)  surf_usm_h%v_0(m) ) * dd2zu(k)2785 2790 dudz(k) = ( u(k+1,j,i)  surf_usm_h%u_0(m) ) * dd2zu(k) 2791 dvdz(k) = ( v(k+1,j,i)  surf_usm_h%v_0(m) ) * dd2zu(k) 2792 2786 2793 ENDDO 2787 2794 ! 2788 ! Compute gradients at downwardfacing walls, only for 2795 ! Compute gradients at downwardfacing walls, only for 2789 2796 ! nonnatural default surfaces 2790 2797 surf_s = surf_def_h(1)%start_index(j,i) … … 2793 2800 k = surf_def_h(1)%k(m) 2794 2801 2795 dudz(k ,j) = ( surf_def_h(1)%u_0(m)  u(k1,j,i) ) * dd2zu(k)2796 dvdz(k ,j) = ( surf_def_h(1)%v_0(m)  v(k1,j,i) ) * dd2zu(k)2802 dudz(k) = ( surf_def_h(1)%u_0(m)  u(k1,j,i) ) * dd2zu(k) 2803 dvdz(k) = ( surf_def_h(1)%v_0(m)  v(k1,j,i) ) * dd2zu(k) 2797 2804 2798 2805 ENDDO 2806 2807 2808 ENDIF 2809 2810 2811 DO k = nzb+1, nzt 2812 2813 def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & 2814 dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & 2815 dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 2816 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + & 2817 dwdy(k)*dvdz(k) ) 2818 2819 IF ( def < 0.0_wp ) def = 0.0_wp 2820 2821 flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) ) 2822 2823 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag 2824 2799 2825 ENDDO 2800 2826 2801 DO j = nys, nyn 2802 DO k = nzb+1, nzt 2803 2804 def = 2.0_wp * ( dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2 ) + & 2805 dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 + & 2806 dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 + & 2807 2.0_wp * ( dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j) + & 2808 dwdy(k,j)*dvdz(k,j) ) 2809 2810 IF ( def < 0.0_wp ) def = 0.0_wp 2811 2812 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 2813 2814 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag 2815 2816 ENDDO 2817 ENDDO 2818 2819 ELSE 2820 2821 DO j = nys, nyn 2822 ! 2823 ! Calculate TKE production by shear. Here, no additional 2824 ! wallbounded code is considered. 2825 ! Why? 2826 DO k = nzb+1, nzt 2827 2828 dudx(k,j) = ( u(k,j,i+1)  u(k,j,i) ) * ddx 2829 dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1)  & 2830 u(k,j1,i)  u(k,j1,i+1) ) * ddy 2831 dudz(k,j) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)  & 2832 u(k1,j,i)  u(k1,j,i+1) ) * & 2833 dd2zu(k) 2834 2835 dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1)  & 2836 v(k,j,i1)  v(k,j+1,i1) ) * ddx 2837 dvdy(k,j) = ( v(k,j+1,i)  v(k,j,i) ) * ddy 2838 dvdz(k,j) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)  & 2839 v(k1,j,i)  v(k1,j+1,i) ) * & 2840 dd2zu(k) 2841 2842 dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k1,j,i+1)  & 2843 w(k,j,i1)  w(k1,j,i1) ) * ddx 2844 dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k1,j+1,i)  & 2845 w(k,j1,i)  w(k1,j1,i) ) * ddy 2846 dwdz(k,j) = ( w(k,j,i)  w(k1,j,i) ) * & 2847 ddzw(k) 2848 2849 def = 2.0_wp * ( & 2850 dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2 & 2851 ) + & 2852 dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 + & 2853 dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 + & 2854 2.0_wp * ( & 2855 dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j) + & 2856 dwdy(k,j)*dvdz(k,j) & 2857 ) 2858 2859 IF ( def < 0.0_wp ) def = 0.0_wp 2860 2861 flag = MERGE( 1.0_wp, 0.0_wp, & 2862 BTEST( wall_flags_0(k,j,i), 29 ) ) 2863 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag 2864 2865 ENDDO 2866 ENDDO 2867 2868 ENDIF 2869 2870 ! 2871 ! If required, calculate TKE production by buoyancy 2872 IF ( .NOT. neutral ) THEN 2873 2874 IF ( .NOT. humidity ) THEN 2875 2876 IF ( ocean_mode ) THEN 2877 ! 2878 ! So far in the ocean no special treatment of density flux 2879 ! in the bottom and top surface layer 2827 2828 ENDDO 2829 ENDDO 2830 2831 ! 2832 ! If required, calculate TKE production by buoyancy 2833 IF ( .NOT. neutral ) THEN 2834 2835 IF ( .NOT. humidity ) THEN 2836 2837 IF ( ocean_mode ) THEN 2838 ! 2839 ! So far in the ocean no special treatment of density flux 2840 ! in the bottom and top surface layer 2841 DO i = nxl, nxr 2880 2842 DO j = nys, nyn 2881 2843 DO k = nzb+1, nzt … … 2891 2853 MERGE( 1.0_wp, 0.0_wp, & 2892 2854 BTEST( wall_flags_0(k,j,i), 9 ) & 2893 ) 2855 ) 2894 2856 ENDDO 2895 2857 ! … … 2921 2883 use_single_reference_value ) * & 2922 2884 drho_air_zw(k) * & 2923 surf_def_h(2)%shf(m) 2885 surf_def_h(2)%shf(m) 2924 2886 ENDDO 2925 2887 ENDIF 2926 2888 2927 2889 ENDDO 2928 2929 ELSE 2930 2890 ENDDO 2891 2892 ELSE ! or IF ( .NOT. ocean_mode ) THEN 2893 2894 DO i = nxl, nxr 2931 2895 DO j = nys, nyn 2896 2932 2897 DO k = nzb+1, nzt 2933 2898 ! … … 2945 2910 MERGE( 1.0_wp, 0.0_wp, & 2946 2911 BTEST( wall_flags_0(k,j,i), 9 ) & 2947 ) 2912 ) 2948 2913 ENDDO 2949 2914 … … 2960 2925 use_single_reference_value ) & 2961 2926 * drho_air_zw(k1) & 2962 * surf_def_h(l)%shf(m) 2963 ENDDO 2927 * surf_def_h(l)%shf(m) 2928 ENDDO 2964 2929 ENDDO 2965 2930 ! … … 2973 2938 use_single_reference_value ) & 2974 2939 * drho_air_zw(k1) & 2975 * surf_lsm_h%shf(m) 2940 * surf_lsm_h%shf(m) 2976 2941 ENDDO 2977 2942 ! … … 2985 2950 use_single_reference_value ) & 2986 2951 * drho_air_zw(k1) & 2987 * surf_usm_h%shf(m) 2988 ENDDO 2952 * surf_usm_h%shf(m) 2953 ENDDO 2989 2954 ENDIF 2990 2955 … … 2998 2963 use_single_reference_value ) & 2999 2964 * drho_air_zw(k) & 3000 * surf_def_h(2)%shf(m) 2965 * surf_def_h(2)%shf(m) 3001 2966 ENDDO 3002 2967 ENDIF 2968 3003 2969 ENDDO 3004 3005 ENDIF 3006 3007 ELSE 3008 2970 ENDDO 2971 2972 ENDIF ! from IF ( .NOT. ocean_mode ) 2973 2974 ELSE ! or IF ( humidity ) THEN 2975 2976 DO i = nxl, nxr 3009 2977 DO j = nys, nyn 3010 2978 … … 3059 3027 k1 = 1.0_wp + 0.61_wp * q(k,j,i)  ql(k,j,i) 3060 3028 k2 = 0.61_wp * pt(k,j,i) 3061 tend(k,j,i) = tend(k,j,i)  & 3029 tend(k,j,i) = tend(k,j,i)  & 3062 3030 kh(k,j,i) * g / & 3063 3031 MERGE( vpt_reference, vpt(k,j,i), & … … 3077 3045 ENDDO 3078 3046 3079 ENDDO 3080 3081 IF ( use_surface_fluxes ) THEN 3082 3083 DO j = nys, nyn 3047 IF ( use_surface_fluxes ) THEN 3084 3048 ! 3085 3049 ! Treat horizontal default surfaces … … 3161 3125 surf_e = surf_usm_h%end_index(j,i) 3162 3126 DO m = surf_s, surf_e 3163 k = surf_ usm_h%k(m)3127 k = surf_lsm_h%k(m) 3164 3128 3165 3129 IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN … … 3193 3157 ENDDO 3194 3158 3195 ENDDO 3196 3197 ENDIF 3198 3199 IF ( use_top_fluxes ) THEN 3200 3201 DO j = nys, nyn 3159 ENDIF ! from IF ( use_surface_fluxes ) THEN 3160 3161 IF ( use_top_fluxes ) THEN 3202 3162 3203 3163 surf_s = surf_def_h(2)%start_index(j,i) … … 3237 3197 ENDDO 3238 3198 3239 ENDDO 3240 3241 ENDIF 3242 3243 ENDIF 3199 ENDIF ! from IF ( use_top_fluxes ) THEN 3200 3201 ENDDO 3202 ENDDO 3244 3203 3245 3204 ENDIF 3246 3205 3247 END DO3206 ENDIF 3248 3207 3249 3208 END SUBROUTINE production_e
Note: See TracChangeset
for help on using the changeset viewer.