Changeset 2118 for palm/trunk/SOURCE/production_e.f90
- Timestamp:
- Jan 17, 2017 4:38:49 PM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/production_e.f90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/production_e.f90
r2101 r2118 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! OpenACC version of subroutine removed 23 23 ! 24 24 ! Former revisions: … … 103 103 104 104 USE wall_fluxes_mod, & 105 ONLY: wall_fluxes_e , wall_fluxes_e_acc105 ONLY: wall_fluxes_e 106 106 107 107 USE kinds 108 108 109 109 PRIVATE 110 PUBLIC production_e, production_e_ acc, production_e_init110 PUBLIC production_e, production_e_init 111 111 112 112 LOGICAL, SAVE :: first_call = .TRUE. !< … … 120 120 END INTERFACE production_e 121 121 122 INTERFACE production_e_acc123 MODULE PROCEDURE production_e_acc124 END INTERFACE production_e_acc125 126 122 INTERFACE production_e_init 127 123 MODULE PROCEDURE production_e_init … … 740 736 ! Description: 741 737 ! ------------ 742 !> Call for all grid points - accelerator version743 !------------------------------------------------------------------------------!744 SUBROUTINE production_e_acc745 746 USE arrays_3d, &747 ONLY: ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho_ocean, shf, &748 tend, tswst, u, v, vpt, w749 750 USE cloud_parameters, &751 ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt752 753 USE control_parameters, &754 ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, &755 humidity, kappa, neutral, ocean, pt_reference, &756 rho_reference, topography, use_single_reference_value, &757 use_surface_fluxes, use_top_fluxes758 759 USE grid_variables, &760 ONLY: ddx, dx, ddy, dy, wall_e_x, wall_e_y761 762 USE indices, &763 ONLY: i_left, i_right, j_north, j_south, nxl, nxr, nys, nyn, nzb, &764 nzb_diff_s_inner, nzb_diff_s_outer, nzb_s_inner, nzt, &765 nzt_diff766 767 IMPLICIT NONE768 769 INTEGER(iwp) :: i !<770 INTEGER(iwp) :: j !<771 INTEGER(iwp) :: k !<772 773 REAL(wp) :: def !<774 REAL(wp) :: dudx !<775 REAL(wp) :: dudy !<776 REAL(wp) :: dudz !<777 REAL(wp) :: dvdx !<778 REAL(wp) :: dvdy !<779 REAL(wp) :: dvdz !<780 REAL(wp) :: dwdx !<781 REAL(wp) :: dwdy !<782 REAL(wp) :: dwdz !<783 REAL(wp) :: k1 !<784 REAL(wp) :: k2 !<785 REAL(wp) :: km_neutral !<786 REAL(wp) :: theta !<787 REAL(wp) :: temp !<788 789 REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs !<790 REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus !<791 REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus !<792 REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsvs !<793 !$acc declare create ( usvs, vsus, wsus, wsvs )794 795 !796 !-- First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at797 !-- vertical walls, if neccessary798 !-- CAUTION: results are slightly different from the ij-version!!799 !-- ij-version should be called further below within the ij-loops!!800 IF ( topography /= 'flat' ) THEN801 CALL wall_fluxes_e_acc( usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, wall_e_y )802 CALL wall_fluxes_e_acc( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, wall_e_y )803 CALL wall_fluxes_e_acc( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, wall_e_x )804 CALL wall_fluxes_e_acc( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, wall_e_x )805 ENDIF806 807 808 !809 !-- Calculate TKE production by shear810 !$acc kernels present( ddzw, dd2zu, kh, km, nzb_diff_s_inner, nzb_diff_s_outer ) &811 !$acc present( nzb_s_inner, pt, q, ql, qsws, qswst, rho_ocean ) &812 !$acc present( shf, tend, tswst, u, v, vpt, w, wall_e_x, wall_e_y ) &813 !$acc copyin( u_0, v_0 )814 DO i = i_left, i_right815 DO j = j_south, j_north816 DO k = 1, nzt817 818 IF ( k >= nzb_diff_s_outer(j,i) ) THEN819 820 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx821 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &822 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy823 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &824 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)825 826 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &827 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx828 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy829 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &830 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)831 832 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &833 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx834 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &835 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy836 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)837 838 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + &839 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &840 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )841 842 IF ( def < 0.0_wp ) def = 0.0_wp843 844 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def845 846 ENDIF847 848 ENDDO849 ENDDO850 ENDDO851 852 IF ( constant_flux_layer ) THEN853 854 !855 !-- Position beneath wall856 !-- (2) - Will allways be executed.857 !-- 'bottom and wall: use u_0,v_0 and wall functions'858 DO i = i_left, i_right859 DO j = j_south, j_north860 DO k = 1, nzt861 862 IF ( ( wall_e_x(j,i) /= 0.0_wp ).OR.( wall_e_y(j,i) /= 0.0_wp ) ) &863 THEN864 865 IF ( k == nzb_diff_s_inner(j,i) - 1 ) THEN866 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx867 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &868 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)869 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy870 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &871 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)872 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)873 874 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN875 !876 !-- Inconsistency removed: as the thermal stratification is877 !-- not taken into account for the evaluation of the wall878 !-- fluxes at vertical walls, the eddy viscosity km must not879 !-- be used for the evaluation of the velocity gradients dudy880 !-- and dwdy881 !-- Note: The validity of the new method has not yet been882 !-- shown, as so far no suitable data for a validation883 !-- has been available884 ! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &885 ! usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )886 ! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &887 ! wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )888 km_neutral = kappa * &889 ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25_wp * &890 0.5_wp * dy891 IF ( km_neutral > 0.0_wp ) THEN892 dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral893 dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral894 ELSE895 dudy = 0.0_wp896 dwdy = 0.0_wp897 ENDIF898 ELSE899 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &900 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy901 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &902 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy903 ENDIF904 905 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN906 !907 !-- Inconsistency removed: as the thermal stratification is908 !-- not taken into account for the evaluation of the wall909 !-- fluxes at vertical walls, the eddy viscosity km must not910 !-- be used for the evaluation of the velocity gradients dvdx911 !-- and dwdx912 !-- Note: The validity of the new method has not yet been913 !-- shown, as so far no suitable data for a validation914 !-- has been available915 ! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &916 ! vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )917 ! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &918 ! wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )919 km_neutral = kappa * &920 ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25_wp * &921 0.5_wp * dx922 IF ( km_neutral > 0.0_wp ) THEN923 dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral924 dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral925 ELSE926 dvdx = 0.0_wp927 dwdx = 0.0_wp928 ENDIF929 ELSE930 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &931 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx932 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &933 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx934 ENDIF935 936 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + &937 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &938 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )939 940 IF ( def < 0.0_wp ) def = 0.0_wp941 942 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def943 944 ENDIF945 !946 !-- (3) - will be executed only, if there is at least one level947 !-- between (2) and (4), i.e. the topography must have a948 !-- minimum height of 2 dz. Wall fluxes for this case have949 !-- already been calculated for (2).950 !-- 'wall only: use wall functions'951 952 IF ( k >= nzb_diff_s_inner(j,i) .AND. &953 k <= nzb_diff_s_outer(j,i)-2 ) THEN954 955 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx956 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &957 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)958 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy959 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &960 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)961 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)962 963 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN964 !965 !-- Inconsistency removed: as the thermal stratification966 !-- is not taken into account for the evaluation of the967 !-- wall fluxes at vertical walls, the eddy viscosity km968 !-- must not be used for the evaluation of the velocity969 !-- gradients dudy and dwdy970 !-- Note: The validity of the new method has not yet971 !-- been shown, as so far no suitable data for a972 !-- validation has been available973 km_neutral = kappa * ( usvs(k,j,i)**2 + &974 wsvs(k,j,i)**2 )**0.25_wp * 0.5_wp * dy975 IF ( km_neutral > 0.0_wp ) THEN976 dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral977 dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral978 ELSE979 dudy = 0.0_wp980 dwdy = 0.0_wp981 ENDIF982 ELSE983 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &984 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy985 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &986 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy987 ENDIF988 989 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN990 !991 !-- Inconsistency removed: as the thermal stratification992 !-- is not taken into account for the evaluation of the993 !-- wall fluxes at vertical walls, the eddy viscosity km994 !-- must not be used for the evaluation of the velocity995 !-- gradients dvdx and dwdx996 !-- Note: The validity of the new method has not yet997 !-- been shown, as so far no suitable data for a998 !-- validation has been available999 km_neutral = kappa * ( vsus(k,j,i)**2 + &1000 wsus(k,j,i)**2 )**0.25_wp * 0.5_wp * dx1001 IF ( km_neutral > 0.0_wp ) THEN1002 dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral1003 dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral1004 ELSE1005 dvdx = 0.0_wp1006 dwdx = 0.0_wp1007 ENDIF1008 ELSE1009 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &1010 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1011 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &1012 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1013 ENDIF1014 1015 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1016 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &1017 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1018 1019 IF ( def < 0.0_wp ) def = 0.0_wp1020 1021 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def1022 1023 ENDIF1024 1025 !1026 !-- (4) - will allways be executed.1027 !-- 'special case: free atmosphere' (as for case (0))1028 IF ( k == nzb_diff_s_outer(j,i)-1 ) THEN1029 1030 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx1031 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &1032 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1033 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &1034 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1035 1036 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &1037 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1038 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1039 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &1040 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1041 1042 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &1043 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1044 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &1045 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1046 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1047 1048 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1049 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &1050 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1051 1052 IF ( def < 0.0_wp ) def = 0.0_wp1053 1054 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def1055 1056 ENDIF1057 1058 ENDIF1059 1060 ENDDO1061 ENDDO1062 ENDDO1063 1064 !1065 !-- Position without adjacent wall1066 !-- (1) - will allways be executed.1067 !-- 'bottom only: use u_0,v_0'1068 DO i = i_left, i_right1069 DO j = j_south, j_north1070 DO k = 1, nzt1071 1072 IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &1073 THEN1074 1075 IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN1076 1077 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx1078 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &1079 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1080 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &1081 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)1082 1083 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &1084 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1085 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1086 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &1087 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)1088 1089 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &1090 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1091 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &1092 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1093 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1094 1095 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1096 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &1097 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1098 1099 IF ( def < 0.0_wp ) def = 0.0_wp1100 1101 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def1102 1103 ENDIF1104 1105 ENDIF1106 1107 ENDDO1108 ENDDO1109 ENDDO1110 1111 ELSEIF ( use_surface_fluxes ) THEN1112 1113 DO i = i_left, i_right1114 DO j = j_south, j_north1115 DO k = 1, nzt1116 1117 IF ( k == nzb_diff_s_outer(j,i)-1 ) THEN1118 1119 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx1120 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &1121 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1122 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &1123 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1124 1125 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &1126 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1127 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1128 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &1129 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1130 1131 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &1132 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1133 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &1134 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1135 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1136 1137 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1138 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &1139 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1140 1141 IF ( def < 0.0_wp ) def = 0.0_wp1142 1143 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def1144 1145 ENDIF1146 1147 ENDDO1148 ENDDO1149 ENDDO1150 1151 ENDIF1152 1153 !1154 !-- If required, calculate TKE production by buoyancy1155 IF ( .NOT. neutral ) THEN1156 1157 IF ( .NOT. humidity ) THEN1158 1159 IF ( use_single_reference_value ) THEN1160 1161 IF ( ocean ) THEN1162 !1163 !-- So far in the ocean no special treatment of density flux1164 !-- in the bottom and top surface layer1165 DO i = i_left, i_right1166 DO j = j_south, j_north1167 DO k = 1, nzt1168 IF ( k > nzb_s_inner(j,i) ) THEN1169 tend(k,j,i) = tend(k,j,i) + &1170 kh(k,j,i) * g / rho_reference * &1171 ( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * &1172 dd2zu(k)1173 ENDIF1174 ENDDO1175 ENDDO1176 ENDDO1177 1178 ELSE1179 1180 DO i = i_left, i_right1181 DO j = j_south, j_north1182 DO k = 1, nzt_diff1183 IF ( k >= nzb_diff_s_inner(j,i) ) THEN1184 tend(k,j,i) = tend(k,j,i) - &1185 kh(k,j,i) * g / pt_reference * &1186 ( pt(k+1,j,i) - pt(k-1,j,i) ) * &1187 dd2zu(k)1188 ENDIF1189 1190 IF ( k == nzb_diff_s_inner(j,i)-1 .AND. &1191 use_surface_fluxes ) THEN1192 tend(k,j,i) = tend(k,j,i) + g / pt_reference * &1193 shf(j,i)1194 ENDIF1195 1196 IF ( k == nzt .AND. use_top_fluxes ) THEN1197 tend(k,j,i) = tend(k,j,i) + g / pt_reference * &1198 tswst(j,i)1199 ENDIF1200 ENDDO1201 ENDDO1202 ENDDO1203 1204 ENDIF1205 1206 ELSE1207 1208 IF ( ocean ) THEN1209 !1210 !-- So far in the ocean no special treatment of density flux1211 !-- in the bottom and top surface layer1212 DO i = i_left, i_right1213 DO j = j_south, j_north1214 DO k = 1, nzt1215 IF ( k > nzb_s_inner(j,i) ) THEN1216 tend(k,j,i) = tend(k,j,i) + &1217 kh(k,j,i) * g / rho_ocean(k,j,i) * &1218 ( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * &1219 dd2zu(k)1220 ENDIF1221 ENDDO1222 ENDDO1223 ENDDO1224 1225 ELSE1226 1227 DO i = i_left, i_right1228 DO j = j_south, j_north1229 DO k = 1, nzt_diff1230 IF( k >= nzb_diff_s_inner(j,i) ) THEN1231 tend(k,j,i) = tend(k,j,i) - &1232 kh(k,j,i) * g / pt(k,j,i) * &1233 ( pt(k+1,j,i) - pt(k-1,j,i) ) * &1234 dd2zu(k)1235 ENDIF1236 1237 IF ( k == nzb_diff_s_inner(j,i)-1 .AND. &1238 use_surface_fluxes ) THEN1239 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &1240 shf(j,i)1241 ENDIF1242 1243 IF ( k == nzt .AND. use_top_fluxes ) THEN1244 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &1245 tswst(j,i)1246 ENDIF1247 ENDDO1248 ENDDO1249 ENDDO1250 1251 ENDIF1252 1253 ENDIF1254 1255 ELSE1256 !1257 !++ This part gives the PGI compiler problems in the previous loop1258 !++ even without any acc statements????1259 ! STOP '+++ production_e problems with acc-directives'1260 ! !acc loop1261 ! DO i = nxl, nxr1262 ! DO j = nys, nyn1263 ! !acc loop vector1264 ! DO k = 1, nzt_diff1265 !1266 ! IF ( k >= nzb_diff_s_inner(j,i) ) THEN1267 !1268 ! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN1269 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i)1270 ! k2 = 0.61_wp * pt(k,j,i)1271 ! tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &1272 ! g / vpt(k,j,i) * &1273 ! ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &1274 ! k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &1275 ! ) * dd2zu(k)1276 ! ELSE IF ( cloud_physics ) THEN1277 ! IF ( ql(k,j,i) == 0.0_wp ) THEN1278 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i)1279 ! k2 = 0.61_wp * pt(k,j,i)1280 ! ELSE1281 ! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)1282 ! temp = theta * t_d_pt(k)1283 ! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &1284 ! ( q(k,j,i) - ql(k,j,i) ) * &1285 ! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &1286 ! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &1287 ! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )1288 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )1289 ! ENDIF1290 ! tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &1291 ! g / vpt(k,j,i) * &1292 ! ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &1293 ! k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &1294 ! ) * dd2zu(k)1295 ! ELSE IF ( cloud_droplets ) THEN1296 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)1297 ! k2 = 0.61_wp * pt(k,j,i)1298 ! tend(k,j,i) = tend(k,j,i) - &1299 ! kh(k,j,i) * g / vpt(k,j,i) * &1300 ! ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + &1301 ! k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &1302 ! pt(k,j,i) * ( ql(k+1,j,i) - &1303 ! ql(k-1,j,i) ) ) * dd2zu(k)1304 ! ENDIF1305 !1306 ! ENDIF1307 !1308 ! ENDDO1309 ! ENDDO1310 ! ENDDO1311 !1312 1313 !!++ Next two loops are probably very inefficiently parallellized1314 !!++ and will require better optimization1315 ! IF ( use_surface_fluxes ) THEN1316 !1317 ! !acc loop1318 ! DO i = nxl, nxr1319 ! DO j = nys, nyn1320 ! !acc loop vector1321 ! DO k = 1, nzt_diff1322 !1323 ! IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN1324 !1325 ! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN1326 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i)1327 ! k2 = 0.61_wp * pt(k,j,i)1328 ! ELSE IF ( cloud_physics ) THEN1329 ! IF ( ql(k,j,i) == 0.0_wp ) THEN1330 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i)1331 ! k2 = 0.61_wp * pt(k,j,i)1332 ! ELSE1333 ! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)1334 ! temp = theta * t_d_pt(k)1335 ! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &1336 ! ( q(k,j,i) - ql(k,j,i) ) * &1337 ! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /&1338 ! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &1339 ! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )1340 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )1341 ! ENDIF1342 ! ELSE IF ( cloud_droplets ) THEN1343 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)1344 ! k2 = 0.61_wp * pt(k,j,i)1345 ! ENDIF1346 !1347 ! tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &1348 ! ( k1* shf(j,i) + k2 * qsws(j,i) )1349 ! ENDIF1350 !1351 ! ENDDO1352 ! ENDDO1353 ! ENDDO1354 !1355 ! ENDIF1356 !1357 ! IF ( use_top_fluxes ) THEN1358 !1359 ! !acc loop1360 ! DO i = nxl, nxr1361 ! DO j = nys, nyn1362 ! !acc loop vector1363 ! DO k = 1, nzt1364 ! IF ( k == nzt ) THEN1365 !1366 ! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN1367 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i)1368 ! k2 = 0.61_wp * pt(k,j,i)1369 ! ELSE IF ( cloud_physics ) THEN1370 ! IF ( ql(k,j,i) == 0.0_wp ) THEN1371 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i)1372 ! k2 = 0.61_wp * pt(k,j,i)1373 ! ELSE1374 ! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)1375 ! temp = theta * t_d_pt(k)1376 ! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &1377 ! ( q(k,j,i) - ql(k,j,i) ) * &1378 ! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /&1379 ! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &1380 ! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )1381 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )1382 ! ENDIF1383 ! ELSE IF ( cloud_droplets ) THEN1384 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)1385 ! k2 = 0.61_wp * pt(k,j,i)1386 ! ENDIF1387 !1388 ! tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &1389 ! ( k1* tswst(j,i) + k2 * qswst(j,i) )1390 !1391 ! ENDIF1392 !1393 ! ENDDO1394 ! ENDDO1395 ! ENDDO1396 !1397 ! ENDIF1398 1399 ENDIF1400 1401 ENDIF1402 !$acc end kernels1403 1404 END SUBROUTINE production_e_acc1405 1406 1407 !------------------------------------------------------------------------------!1408 ! Description:1409 ! ------------1410 738 !> Call for grid point i,j 1411 739 !------------------------------------------------------------------------------!
Note: See TracChangeset
for help on using the changeset viewer.
