Changeset 1746 for trunk/MESOSCALE/LMD_MM_MARS/SRC
- Timestamp:
- Jul 24, 2017, 5:06:29 PM (7 years ago)
- Location:
- trunk/MESOSCALE/LMD_MM_MARS/SRC
- Files:
-
- 2 edited
-
LES/modif_mars/module_initialize_les.F (modified) (24 diffs)
-
WRFV2/dyn_em/module_initialize_real.F (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/MESOSCALE/LMD_MM_MARS/SRC/LES/modif_mars/module_initialize_les.F
r1724 r1746 86 86 87 87 INTEGER, PARAMETER :: nl_max = 1000 88 REAL, DIMENSION(nl_max) :: zk, p_in, theta, Tz,rho, u, v, qv, pd_in88 REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in 89 89 INTEGER :: nl_in 90 90 … … 116 116 ! for mode 3 117 117 REAL, DIMENSION(nl_max) :: profdustq,profdustn 118 REAL, DIMENSION(nl_max) :: hrsw,hrlw,hrdyn 119 REAL, DIMENSION(nl_max) :: trac 120 REAL, DIMENSION(nl_max) :: prescribed_sw,prescribed_lw,prescribed_dyn 121 REAL, DIMENSION(nl_max) :: venus_hrdyn 122 REAL, DIMENSION(nl_max) :: altitude 118 REAL, DIMENSION(nl_max) :: prescribed_sw,prescribed_lw 123 119 !!MARS 124 120 … … 127 123 !INTEGER :: hypsometric_opt = 2 ! Wee et al. 2012 correction 128 124 125 character(len=10) :: planet 129 126 130 127 #ifdef DM_PARALLEL 131 128 # include <data_calls.inc> 132 129 #endif 130 131 call init_planet_constants 133 132 134 133 SELECT CASE ( model_data_order ) … … 173 172 174 173 175 stretch_grid = .false.174 !stretch_grid = .false. 176 175 ! FOR LES, set stretch to false 177 !stretch_grid = .true.176 stretch_grid = .true. 178 177 delt = 3. 179 178 ! z_scale = .50 180 179 ! z_scale = .10 181 180 ! z_scale = .25 182 ! z_scale = .1 181 ! z_scale = .15 183 182 pi = 2.*asin(1.0) 184 183 write(6,*) ' pi is ',pi … … 303 302 grid%rdx = 1./config_flags%dx 304 303 grid%rdy = 1./config_flags%dy 304 305 305 ! get the sounding from the ascii sounding file, first get dry sounding and 306 306 ! calculate base state … … 310 310 write(6,*) ' getting dry sounding for base state ' 311 311 312 CALL get_sounding( zk, p_in, pd_in, theta, Tz,rho, u, v, qv, dry_sounding, nl_max, nl_in )312 CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in ) 313 313 ENDIF 314 314 CALL wrf_dm_bcast_real( zk , nl_max ) … … 316 316 CALL wrf_dm_bcast_real( pd_in , nl_max ) 317 317 CALL wrf_dm_bcast_real( theta , nl_max ) 318 CALL wrf_dm_bcast_real( Tz , nl_max )319 318 CALL wrf_dm_bcast_real( rho , nl_max ) 320 319 CALL wrf_dm_bcast_real( u , nl_max ) … … 385 384 ENDDO 386 385 387 xs=ide/2 -25 388 !xs=ids -3 389 !xe=xs + 6 390 xe=xs + 25 391 !ys=jde/2 -3 392 ys=jde/2 -25 393 !ye=ys + 6 394 ye=ys + 25 395 mtn_ht = 1000 386 xs=ide/2 -3 387 xs=ids -3 388 xe=xs + 6 389 ys=jde/2 -3 390 ye=ys + 6 391 mtn_ht = 500 396 392 #ifdef MTN 397 393 DO j=max(ys,jds),min(ye,jde-1) 398 394 DO i=max(xs,ids),min(xe,ide-1) 399 grid%ht(i,j) = alt_input +mtn_ht * 0.25 * &395 grid%ht(i,j) = mtn_ht * 0.25 * & 400 396 ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) ) * & 401 397 ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) ) … … 438 434 p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top 439 435 grid%pb(i,k,j) = p_level 440 grid%t_init(i,k,j) = interp_0_log( theta, p_in, p_level, nl_in ) - t0 441 grid%tz_init(i,k,j) = interp_0_log( Tz, p_in, p_level, nl_in ) 436 grid%t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0 442 437 grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm 443 438 ENDDO … … 463 458 ENDDO 464 459 ENDDO 460 465 461 IF ( wrf_dm_on_monitor() ) THEN 466 462 write(6,*) ' ptop is ',grid%p_top … … 472 468 write(6,*) ' getting moist sounding for full state ' 473 469 dry_sounding = .false. 474 CALL get_sounding( zk, p_in, pd_in, theta, Tz,rho, u, v, qv, dry_sounding, nl_max, nl_in )470 CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in ) 475 471 476 472 DO J = jts, min(jde-1,jte) … … 495 491 p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top 496 492 497 moist(i,k,j,P_QV) = interp_0 _log( qv, pd_in, p_level, nl_in )498 grid%t_1(i,k,j) = interp_0 _log( theta, pd_in, p_level, nl_in ) - t0493 moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in ) 494 grid%t_1(i,k,j) = interp_0( theta, pd_in, p_level, nl_in ) - t0 499 495 grid%t_2(i,k,j) = grid%t_1(i,k,j) 500 496 … … 698 694 DO K = 1, kte-1 699 695 p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top 700 grid%v_1(i,k,j) = interp_0 _log( v, p_in, p_level, nl_in )696 grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in ) 701 697 grid%v_2(i,k,j) = grid%v_1(i,k,j) 702 698 ENDDO … … 722 718 DO K = 1, kte-1 723 719 p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top 724 grid%u_1(i,k,j) = interp_0 _log( u, p_in, p_level, nl_in )720 grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in ) 725 721 grid%u_2(i,k,j) = grid%u_1(i,k,j) 726 722 ENDDO … … 851 847 DO k=1,kte-1 852 848 p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top 853 scalar(its:ite,k,jts:jte,2) = interp_0 _log( qv, pd_in, p_level, nl_in )849 scalar(its:ite,k,jts:jte,2) = interp_0( qv, pd_in, p_level, nl_in ) 854 850 scalar(its:ite,k,jts:jte,3) = 0. 855 851 !! water ice is set to 0 (was put into water vapor when building prof from MCD) … … 866 862 DO k=1,kte-1 867 863 p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top 868 scalar(its:ite,k,jts:jte,4) = interp_0 _log( profdustq, pd_in, p_level, nl_in )869 scalar(its:ite,k,jts:jte,5) = interp_0 _log( profdustn, pd_in, p_level, nl_in )864 scalar(its:ite,k,jts:jte,4) = interp_0( profdustq, pd_in, p_level, nl_in ) 865 scalar(its:ite,k,jts:jte,5) = interp_0( profdustn, pd_in, p_level, nl_in ) 870 866 ENDDO 871 867 print *, "DUST Q", scalar(its,:,jts,4) … … 886 882 DO i = its, ite 887 883 !!! we use Q2 as a vehicle for heating rates! sick! 888 grid%m_q2(i,k,j) = interp_0 _log( profdustq, pd_in, p_level, nl_in ) &889 + interp_0 _log( profdustn, pd_in, p_level, nl_in )884 grid%m_q2(i,k,j) = interp_0( profdustq, pd_in, p_level, nl_in ) & 885 + interp_0( profdustn, pd_in, p_level, nl_in ) 890 886 ENDDO 891 887 ENDDO … … 893 889 !print*,k,grid%m_q2(1,k,1) 894 890 ENDDO 895 891 print*,'planet',planet 896 892 IF (planet.eq."prescribed") Then 897 call read_hr( hrsw,hrlw,hrdyn,nl_in)893 call read_hr(profdustq,profdustn,nl_in) 898 894 open(unit=17,file="prescribed_sw.txt",action="write") 899 895 open(unit=18,file="prescribed_lw.txt",action="write") 900 open(unit=19,file="prescribed_dyn.txt",action="write")901 896 DO k=1,kte!-1 902 897 p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top 903 prescribed_sw(k) = interp_0_log( hrsw, pd_in, p_level, nl_in ) 904 prescribed_lw(k) = interp_0_log( hrlw, pd_in, p_level, nl_in ) 905 prescribed_dyn(k) = interp_0_log( hrdyn, pd_in, p_level, nl_in ) 898 prescribed_sw(k) = interp_0( profdustq, pd_in, p_level, nl_in ) 899 prescribed_lw(k) = interp_0( profdustn, pd_in, p_level, nl_in ) 906 900 write (17,*) prescribed_sw(k) 907 901 write (18,*) prescribed_lw(k) 908 write (19,*) prescribed_dyn(k)909 902 ENDDO 910 close(unit=1 9)903 close(unit=17) 911 904 close(unit=18) 912 close(unit=17)913 905 ENDIF 914 915 IF (planet.eq."venus") Then916 call read_hr(hrsw,hrlw,hrdyn,nl_in)917 open(unit=20,file="venus_hrdyn.txt",action="write")918 DO k=1,kte!-1919 p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top920 venus_hrdyn(k) = interp_0_log( hrdyn, pd_in, p_level, nl_in )921 write (20,*) venus_hrdyn(k)922 ENDDO923 close(unit=20)924 ENDIF925 926 open(unit=21,file="altitude.txt",action="write")927 DO k=1,kte!-1928 p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top929 altitude(k) = interp_0_log( zk, pd_in, p_level, nl_in )930 write (21,*) altitude(k)931 ENDDO932 close(unit=21)933 906 934 907 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 968 941 !--------------------------------------------------------------------------- 969 942 970 subroutine get_sounding( zk, p, p_dry, theta, Tz,rho, &943 subroutine get_sounding( zk, p, p_dry, theta, rho, & 971 944 u, v, qv, dry, nl_max, nl_in ) 972 945 implicit none … … 974 947 integer nl_max, nl_in 975 948 real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), & 976 u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max) ,Tz(nl_max)949 u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max) 977 950 logical dry 978 951 … … 1139 1112 v(k) = v_input(k) 1140 1113 qv(k) = qv_input(k) 1141 Tz(k) = t_therm(k)1142 1114 1143 1115 enddo … … 1226 1198 implicit none 1227 1199 integer n 1228 real pdustq(n +1),pdustn(n+1)1200 real pdustq(n),pdustn(n) 1229 1201 logical end_of_file 1230 1202 … … 1251 1223 end subroutine read_dust 1252 1224 !!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1253 subroutine read_hr(hr_sw,hr_lw, hr_dyn,n)1225 subroutine read_hr(hr_sw,hr_lw,n) 1254 1226 implicit none 1255 1227 integer n 1256 real hr_sw(n +1),hr_lw(n+1),hr_dyn(n+1)1228 real hr_sw(n),hr_lw(n) 1257 1229 logical end_of_file 1258 1230 … … 1267 1239 do while (.not. end_of_file) 1268 1240 1269 read(11,*,end=103) hr_sw(k+1),hr_lw(k+1) ,hr_dyn(k+1)1270 write(*,*) k,hr_sw(k+1),hr_lw(k+1) ,hr_dyn(k+1)1241 read(11,*,end=103) hr_sw(k+1),hr_lw(k+1) 1242 write(*,*) k,hr_sw(k+1),hr_lw(k+1) 1271 1243 k = k+1 1272 go to 11 41244 go to 113 1273 1245 103 end_of_file = .true. 1274 1246 114 continue -
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_real.F
r1724 r1746 21 21 USE module_soil_pre 22 22 USE module_date_time 23 use module_init_utilities24 23 #ifdef DM_PARALLEL 25 24 USE module_dm … … 119 118 REAL :: max_dz 120 119 121 INTEGER , PARAMETER :: nl_max = 1000 122 INTEGER :: nl_in 120 ! INTEGER , PARAMETER :: nl_max = 1000 123 121 ! REAL , DIMENSION(nl_max) :: grid%em_dn 124 REAL, DIMENSION(nl_max) :: pd_in 125 REAL, DIMENSION(nl_max) :: profdustq,profdustn 122 126 123 integer::oops1,oops2 127 124 … … 152 149 !LOGICAL :: interp_theta = .false. ! Wee et al. 2012 correction 153 150 REAL :: pfu, pfd, phm 154 REAL :: tpot 151 155 152 156 153 #ifdef DM_PARALLEL … … 1083 1080 !!END DO 1084 1081 1085 IF ( planet == "mars" ) then 1082 1086 1083 !--get vertical size of the GCM input array and allocate new stuff 1087 1088 1089 1090 1091 1092 1093 1084 sizegcm=SIZE(grid%em_rh_gc(1,:,1)) 1085 ALLOCATE(sig(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1)) 1086 !ALLOCATE(ap(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1)) 1087 ALLOCATE(bp(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1)) 1088 1089 DO j = jts , MIN ( jde-1 , jte ) 1090 DO i = its , MIN (ide-1 , ite ) 1094 1091 1095 1092 !!! Define old sigma levels for each column 1096 1093 sig(i,:,j)=grid%em_p_gc(i,:,j)/grid%em_psfc_gc(i,j) 1097 1094 1098 1095 !!! Compute new sigma levels from old sigma levels with GCM (low-res) and WRF (hi-res) surface pressure 1099 1096 !!! (dimlevs,sigma_gcm, ps_gcm, ps_hr, sigma_hr) 1100 1097 CALL build_sigma_hr(sizegcm,sig(i,:,j),grid%em_psfc_gc(i,j),grid%psfc(i,j),bp(i,:,j)) 1101 1098 1102 1099 !!! Calculate new pressure levels 1103 1104 1105 1106 1107 1108 1109 1100 grid%em_pd_gc(i,:,j)=bp(i,:,j)*grid%psfc(i,j) 1101 1102 END DO 1103 END DO 1104 1105 DEALLOCATE(sig) 1106 DEALLOCATE(bp) 1110 1107 1111 1108 !!****MARS who knows... 1112 1109 grid%em_rh_gc(:,:,:)=0. 1113 1110 1114 1111 … … 1116 1113 !grid%em_pd_gc=grid%em_p_gc 1117 1114 !!****MARS 1118 ELSE ! VENUS 1115 1119 1116 1120 1117 … … 1123 1120 !! dry top pressure (constant). 1124 1121 ! 1125 CALL p_dts ( grid%em_mu0 , grid%em_intq_gc , grid%psfc , grid%p_top , & 1126 ids , ide , jds , jde , 1 , num_metgrid_levels , & 1127 ims , ime , jms , jme , 1 , num_metgrid_levels , & 1128 its , ite , jts , jte , 1 , num_metgrid_levels ) 1129 ENDIF 1130 IF ( planet == "mars" ) then 1122 !CALL p_dts ( grid%em_mu0 , grid%em_intq_gc , grid%psfc , grid%p_top , & 1123 ! ids , ide , jds , jde , 1 , num_metgrid_levels , & 1124 ! ims , ime , jms , jme , 1 , num_metgrid_levels , & 1125 ! its , ite , jts , jte , 1 , num_metgrid_levels ) 1126 1131 1127 !!****MARS 1132 1133 1134 1135 1136 1137 1138 1128 DO j = jts , MIN ( jde-1 , jte ) 1129 DO i = its , MIN (ide-1 , ite ) 1130 1131 grid%em_mu0(i,j) = grid%psfc(i,j) - grid%p_top 1132 1133 END DO 1134 END DO 1139 1135 !!****MARS 1140 ELSE ! VENUS 1136 1141 1137 1142 1138 !! Compute the dry, hydrostatic surface pressure. 1143 1139 ! 1144 CALL p_dhs ( grid%em_pdhs , grid%ht , p00 , t00 , a , & 1145 ids , ide , jds , jde , kds , kde , & 1146 ims , ime , jms , jme , kms , kme , & 1147 its , ite , jts , jte , kts , kte ) 1148 ENDIF 1140 !CALL p_dhs ( grid%em_pdhs , grid%ht , p00 , t00 , a , & 1141 ! ids , ide , jds , jde , kds , kde , & 1142 ! ims , ime , jms , jme , kms , kme , & 1143 ! its , ite , jts , jte , kts , kte ) 1149 1144 !!****MARS: voir remarques dans la routine 1150 1145 !!****MARS: dry hydrostatic pressure comes from the GCM ... … … 2299 2294 IF (( i .EQ. its ) .AND. ( j .EQ. jts )) print *, temp, k 2300 2295 !!! MODIF WRFV3.1 - parameter tiso 2301 IF (planet .eq. "mars" ) THEN 2302 grid%em_t_init(i,k,j) = temp*(p00/grid%em_pb(i,k,j))**(r_d/cp) - t0 2303 ELSE 2304 grid%em_t_init(i,k,j) = (temp**nu + nu*(TT00**nu)*log((p00/grid%em_pb(i,k,j))**rcp))**(1/nu) -t0 2305 ENDIF 2296 grid%em_t_init(i,k,j) = temp*(p00/grid%em_pb(i,k,j))**(r_d/cp) - t0 2306 2297 grid%em_alb(i,k,j) = (r_d/p1000mb)*(grid%em_t_init(i,k,j)+t0)*(grid%em_pb(i,k,j)/p1000mb)**cvpm 2307 2298 END DO … … 2610 2601 !!NB: q2 is used for other purpose ... 2611 2602 2603 2612 2604 ! END IF 2613 2605 … … 2788 2780 temp1 = MAX(tiso,t00+A*LOG(grid%em_pb(i,k,j)/p00)) 2789 2781 temp2 = MAX(tiso,t00+A*LOG( pb_int/p00)) 2790 IF (planet .eq. "mars" ) THEN 2791 grid%em_t_init(i,k,j) = temp1*(p00/grid%em_pb(i,k,j))**(r_d/cp) - t0 2792 t_init_int(i,k,j) = temp2*(p00/pb_int )**(r_d/cp) - t0 2793 ELSE 2794 ENDIF 2782 grid%em_t_init(i,k,j) = temp1*(p00/grid%em_pb(i,k,j))**(r_d/cp) - t0 2783 t_init_int(i,k,j) = temp2*(p00/pb_int )**(r_d/cp) - t0 2795 2784 grid%em_alb(i,k,j) = (r_d/p1000mb)*(grid%em_t_init(i,k,j)+t0)*(grid%em_pb(i,k,j)/p1000mb)**cvpm 2796 2785 END DO … … 4980 4969 ! temp = t00 + A*LOG(pb/p00) 4981 4970 temp = MAX ( tiso, t00 + A*LOG(pb/p00) ) 4982 IF (planet .eq. "mars" ) THEN 4983 t_init = temp*(p00/pb)**(r_d/cp) - t0 4984 ELSE 4985 t_init = (temp**nu + nu*(TT00**nu)*log((p00/pb)**(rcp)))**(1/nu) - t0 4986 ENDIF 4971 t_init = temp*(p00/pb)**(r_d/cp) - t0 4987 4972 alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm 4988 4973 END DO … … 5020 5005 ! temp = t00 + A*LOG(pb/p00) 5021 5006 temp = MAX ( tiso, t00 + A*LOG(pb/p00) ) 5022 IF (planet .eq. "mars" ) THEN 5023 t_init = temp*(p00/pb)**(r_d/cp) - t0 5024 ELSE 5025 t_init = (temp**nu + nu*(TT00**nu)*log((p00/pb)**(rcp)))**(1/nu) -t0 5026 ENDIF 5007 t_init = temp*(p00/pb)**(r_d/cp) - t0 5027 5008 alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm 5028 5009 znw(k+1) = znw(k) - dz*g / ( mub*alb(k) ) … … 5044 5025 ! temp = t00 + A*LOG(pb/p00) 5045 5026 temp = MAX ( tiso, t00 + A*LOG(pb/p00) ) 5046 IF (planet .eq. "mars" ) THEN 5047 t_init = temp*(p00/pb)**(r_d/cp) - t0 5048 ELSE 5049 t_init = (temp**nu + nu*(TT00**nu)*log((p00/pb)**(rcp)))**(1/nu) -t0 5050 ENDIF 5027 t_init = temp*(p00/pb)**(r_d/cp) - t0 5051 5028 alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm 5052 5029 znw(k+1) = znw(k) - dz*g / ( mub*alb(k) ) … … 5075 5052 ! temp = t00 + A*LOG(pb/p00) 5076 5053 temp = MAX ( tiso, t00 + A*LOG(pb/p00) ) 5077 IF (planet .eq. "mars" ) THEN 5078 t_init = temp*(p00/pb)**(r_d/cp) - t0 5079 ELSE 5080 t_init = (temp**nu + nu*(TT00**nu)*log((p00/pb)**(rcp)))**(1/nu) -t0 5081 ENDIF 5054 t_init = temp*(p00/pb)**(r_d/cp) - t0 5082 5055 alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm 5083 5056 END DO … … 5747 5720 end subroutine build_sigma_hr 5748 5721 5722 5723 5724 5725 5749 5726 END MODULE module_initialize 5750 5727
Note: See TracChangeset
for help on using the changeset viewer.
