Changeset 1724 for trunk/MESOSCALE/LMD_MM_MARS/SRC/LES

Timestamp:

Jul 21, 2017, 4:07:38 PM (8 years ago)

Author:

mlefevre

Message:

Add of Venus mesoscale structure + update of Venus LES

Location:

trunk/MESOSCALE/LMD_MM_MARS/SRC/LES/modif_mars

Files:

• Registry.EM (modified) (1 diff)
• Registry.EM.newphys (modified) (3 diffs)
• module_big_step_utilities_em.F (modified) (2 diffs)
• module_first_rk_step_part1.F (modified) (1 diff)
• module_initialize_les.F (modified) (24 diffs)

trunk/MESOSCALE/LMD_MM_MARS/SRC/LES/modif_mars/Registry.EM

@@ -99 +99 @@
 state  real  RICE      ikj  misc  1  -  rd   "RICE"      "ICE RADIUS"                      "m"       #SAVEMARS3 rice
 state  real  PDTZ      ikj  misc  1  -  rd   "PDT"       "TEMP TENDENCY"                   "K s-1"   #SAVEMARS3 pdt
+state  real  HR_SW     ikj   misc  1  -  -  "HR_SW"     "HEATING RATE SW"                 "K/s"
+state  real  HR_LW     ikj   misc  1  -  -  "HR_LW"     "HEATING RATE LW"                 "K/s"
+state  real  HR_DYN     ikj   misc  1  -  -  "HR_DYN"     "HEATING RATE DYN"                 "K/s"
+state  real  DDT     ikj   misc  1  -  -  "DDT"     "TEMPERATURE TENDENCY"                 "K/s"
+state  real  DT_RAD     ikj   misc  1  -  -  "DT_RAD"     "RADIATIVE TEMPERATURE TENDENCY"                 "K/s"
+state  real  DT_VDF     ikj   misc  1  -  -  "DT_VDF"     "PBL CLMAIN TEMPERATURE TENDENCY"                 "K/s"
+state  real  DT_AJS     ikj   misc  1  -  -  "DT_AJS"     "AJUSTEMENT SEC TEMPERATURE TENDENCY"                 "K/s"
 ####
 ####

trunk/MESOSCALE/LMD_MM_MARS/SRC/LES/modif_mars/Registry.EM.newphys

@@ -82 +82 @@
 state  real  HR_SW     ikj   misc  1  -  -  "HR_SW"     "HEATING RATE SW"                 "K/s"
 state  real  HR_LW     ikj   misc  1  -  -  "HR_LW"     "HEATING RATE LW"                 "K/s"
+state  real  HR_DYN     ikj   misc  1  -  -  "HR_DYN"     "HEATING RATE DYN"                 "K/s"
+state  real  DDT     ikj   misc  1  -  -  "DDT"     "TEMPERATURE TENDENCY"                 "K/s"
+state  real  DT_RAD     ikj   misc  1  -  -  "DT_RAD"     "RADIATIVE TEMPERATURE TENDENCY"                 "K/s"
+state  real  DT_VDF     ikj   misc  1  -  -  "DT_VDF"     "PBL CLMAIN TEMPERATURE TENDENCY"                 "K/s"
+state  real  DT_AJS     ikj   misc  1  -  -  "DT_AJS"     "AJUSTEMENT SEC TEMPERATURE TENDENCY"                 "K/s"
+
 state  real  SWDOWNZ    ij   misc  1  -  -  "SWDOWNZ"   "DOWNWARD SW FLUX AT SURFACE"     "W m-2"
 #### DUST CYCLE
@@ -96 +102 @@
 state  real  ZMAX       ij   misc  1  -  -  "ZMAX"      "MAXIMUM LEVEL THERMALS"          "m"
 state  real  HFMAX      ij   misc  1  -  -  "HFMAX"     "MAXIMUM PBL HEAT FLUX"           "m.K/s"
+
 ####
 ####
@@ -303 +310 @@
 state    real   t_save         ikj     dyn_em      1         -               "t_save"
 #                                               
-
+state    real   tz_init        ikj     dyn_em      1         -      ir       "tz_init" "initial potential temperature" "K"
 
 # Mass

trunk/MESOSCALE/LMD_MM_MARS/SRC/LES/modif_mars/module_big_step_utilities_em.F

@@ -14 +14 @@
 
    USE module_domain, ONLY : domain
+   USE module_configure
    USE module_model_constants
    USE module_state_description
@@ -4475 +4476 @@
       p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
       pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
-      t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
+      IF (planet.eq. "venus" ) THEN 
+      !call tpot2t(1,th_phy(i,k,j), t_phy(i,k,j), pi_phy(i,k,j))
+        t_phy(i,k,j)= (th_phy(i,k,j)**nu - nu*(TT00**nu)*log((p1000mb/p_phy(i,k,j))**rcp))**(1/nu)
+      ELSE
+        t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
+      ENDIF
       rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
       mu_3d(i,k,j) = mu(i,j)

trunk/MESOSCALE/LMD_MM_MARS/SRC/LES/modif_mars/module_first_rk_step_part1.F

@@ -650 +650 @@
      &       ,HR_SW=grid%HR_SW   &
      &       ,HR_LW=grid%HR_LW   &
+     &       ,HR_DYN=grid%HR_DYN   &
+     &       ,DDT=grid%DDT  &
+     &       ,DT_RAD=grid%DT_RAD  &
+     &       ,DT_VDF=grid%DT_VDF  &
+     &       ,DT_AJS=grid%DT_AJS  &
      &       ,SWDOWNZ=grid%SWDOWNZ   &
      &       ,TAU_DUST=grid%TAU_DUST   &

trunk/MESOSCALE/LMD_MM_MARS/SRC/LES/modif_mars/module_initialize_les.F

@@ -86 +86 @@
 
    INTEGER, PARAMETER :: nl_max = 1000
-   REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in
+   REAL, DIMENSION(nl_max) :: zk, p_in, theta, Tz, rho, u, v, qv, pd_in
    INTEGER :: nl_in
 
@@ -116 +116 @@
  ! for mode 3
  REAL, DIMENSION(nl_max) :: profdustq,profdustn
- REAL, DIMENSION(nl_max) :: prescribed_sw,prescribed_lw
+ REAL, DIMENSION(nl_max) :: hrsw,hrlw,hrdyn
+ REAL, DIMENSION(nl_max) :: trac
+ REAL, DIMENSION(nl_max) :: prescribed_sw,prescribed_lw,prescribed_dyn
+ REAL, DIMENSION(nl_max) :: venus_hrdyn
+ REAL, DIMENSION(nl_max) :: altitude
 !!MARS
 
@@ -123 +127 @@
       !INTEGER :: hypsometric_opt = 2 ! Wee et al. 2012 correction
 
-character(len=10) :: planet
 
 #ifdef DM_PARALLEL
 #    include <data_calls.inc>
 #endif
- 
-   call init_planet_constants
 
    SELECT CASE ( model_data_order )
@@ -172 +173 @@
 
 
-   !stretch_grid = .false.
+   stretch_grid = .false.
 !  FOR LES, set stretch to false
-   stretch_grid = .true.
+   !stretch_grid = .true.
    delt = 3.
 !   z_scale = .50
 !   z_scale = .10
 !   z_scale = .25
-!   z_scale = .15
+!   z_scale = .1
    pi = 2.*asin(1.0)
    write(6,*) ' pi is ',pi
@@ -302 +303 @@
    grid%rdx = 1./config_flags%dx
    grid%rdy = 1./config_flags%dy
-
 !  get the sounding from the ascii sounding file, first get dry sounding and 
 !  calculate base state
@@ -310 +310 @@
   write(6,*) ' getting dry sounding for base state '
 
-  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
+  CALL get_sounding( zk, p_in, pd_in, theta, Tz, rho, u, v, qv, dry_sounding, nl_max, nl_in )
   ENDIF
   CALL wrf_dm_bcast_real( zk , nl_max )
@@ -316 +316 @@
   CALL wrf_dm_bcast_real( pd_in , nl_max )
   CALL wrf_dm_bcast_real( theta , nl_max )
+  CALL wrf_dm_bcast_real( Tz , nl_max )
   CALL wrf_dm_bcast_real( rho , nl_max )
   CALL wrf_dm_bcast_real( u , nl_max )
@@ -384 +385 @@
   ENDDO
 
-  xs=ide/2 -3
-  xs=ids   -3
-  xe=xs + 6
-  ys=jde/2 -3
-  ye=ys + 6
-  mtn_ht = 500
+  xs=ide/2 -25
+  !xs=ids   -3
+  !xe=xs + 6
+  xe=xs + 25
+  !ys=jde/2 -3
+  ys=jde/2 -25
+  !ye=ys + 6
+  ye=ys + 25
+  mtn_ht = 1000
 #ifdef MTN
   DO j=max(ys,jds),min(ye,jde-1)
   DO i=max(xs,ids),min(xe,ide-1)
-     grid%ht(i,j) = mtn_ht * 0.25 * &
+     grid%ht(i,j) = alt_input + mtn_ht * 0.25 * &
                ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) ) * &
                ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) )
@@ -434 +438 @@
       p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
       grid%pb(i,k,j) = p_level
-      grid%t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0
+      grid%t_init(i,k,j) = interp_0_log( theta, p_in, p_level, nl_in ) - t0
+      grid%tz_init(i,k,j) = interp_0_log( Tz, p_in, p_level, nl_in )
       grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
     ENDDO
@@ -458 +463 @@
   ENDDO
   ENDDO
-
   IF ( wrf_dm_on_monitor() ) THEN
     write(6,*) ' ptop is ',grid%p_top
@@ -468 +472 @@
   write(6,*) ' getting moist sounding for full state '
   dry_sounding = .false.
-  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
+  CALL get_sounding( zk, p_in, pd_in, theta, Tz, rho, u, v, qv, dry_sounding, nl_max, nl_in )
 
   DO J = jts, min(jde-1,jte)
@@ -491 +495 @@
       p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top
 
-      moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
-      grid%t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
+      moist(i,k,j,P_QV) = interp_0_log( qv, pd_in, p_level, nl_in )
+      grid%t_1(i,k,j)          = interp_0_log( theta, pd_in, p_level, nl_in ) - t0
       grid%t_2(i,k,j)          = grid%t_1(i,k,j)
       
@@ -694 +698 @@
     DO K = 1, kte-1
       p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
-      grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
+      grid%v_1(i,k,j) = interp_0_log( v, p_in, p_level, nl_in )
       grid%v_2(i,k,j) = grid%v_1(i,k,j)
     ENDDO
@@ -718 +722 @@
     DO K = 1, kte-1
       p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
-      grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
+      grid%u_1(i,k,j) = interp_0_log( u, p_in, p_level, nl_in )
       grid%u_2(i,k,j) = grid%u_1(i,k,j)
     ENDDO
@@ -847 +851 @@
       DO k=1,kte-1
          p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top
-         scalar(its:ite,k,jts:jte,2) = interp_0( qv, pd_in, p_level, nl_in )
+         scalar(its:ite,k,jts:jte,2) = interp_0_log( qv, pd_in, p_level, nl_in )
          scalar(its:ite,k,jts:jte,3) = 0.
            !! water ice is set to 0 (was put into water vapor when building prof from MCD)
@@ -862 +866 @@
       DO k=1,kte-1
          p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top
-         scalar(its:ite,k,jts:jte,4) = interp_0( profdustq, pd_in, p_level, nl_in )
-         scalar(its:ite,k,jts:jte,5) = interp_0( profdustn, pd_in, p_level, nl_in )
+         scalar(its:ite,k,jts:jte,4) = interp_0_log( profdustq, pd_in, p_level, nl_in )
+         scalar(its:ite,k,jts:jte,5) = interp_0_log( profdustn, pd_in, p_level, nl_in )
       ENDDO
       print *, "DUST Q", scalar(its,:,jts,4)
@@ -882 +886 @@
        DO i = its, ite
          !!! we use Q2 as a vehicle for heating rates! sick!
-         grid%m_q2(i,k,j) = interp_0( profdustq, pd_in, p_level, nl_in ) &
-                          + interp_0( profdustn, pd_in, p_level, nl_in )
+         grid%m_q2(i,k,j) = interp_0_log( profdustq, pd_in, p_level, nl_in ) &
+                          + interp_0_log( profdustn, pd_in, p_level, nl_in )
        ENDDO
        ENDDO
@@ -889 +893 @@
        !print*,k,grid%m_q2(1,k,1)
     ENDDO
-    print*,'planet',planet   
+
     IF (planet.eq."prescribed") Then
-      call read_hr(profdustq,profdustn,nl_in)
+      call read_hr(hrsw,hrlw,hrdyn,nl_in)
       open(unit=17,file="prescribed_sw.txt",action="write")
       open(unit=18,file="prescribed_lw.txt",action="write")
+      open(unit=19,file="prescribed_dyn.txt",action="write")
       DO k=1,kte!-1
         p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top
-        prescribed_sw(k) = interp_0( profdustq, pd_in, p_level, nl_in )
-        prescribed_lw(k) = interp_0( profdustn, pd_in, p_level, nl_in )
+        prescribed_sw(k) = interp_0_log( hrsw, pd_in, p_level, nl_in )
+        prescribed_lw(k) = interp_0_log( hrlw, pd_in, p_level, nl_in )
+        prescribed_dyn(k) = interp_0_log( hrdyn, pd_in, p_level, nl_in )
         write (17,*) prescribed_sw(k)
         write (18,*) prescribed_lw(k)
+        write (19,*) prescribed_dyn(k)
       ENDDO
+      close(unit=19)
+      close(unit=18)
       close(unit=17)
-      close(unit=18)
     ENDIF
+    
+    IF (planet.eq."venus") Then
+      call read_hr(hrsw,hrlw,hrdyn,nl_in)
+      open(unit=20,file="venus_hrdyn.txt",action="write")
+      DO k=1,kte!-1
+        p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top
+        venus_hrdyn(k) = interp_0_log( hrdyn, pd_in, p_level, nl_in )
+        write (20,*) venus_hrdyn(k)
+      ENDDO
+      close(unit=20)
+    ENDIF
+
+    open(unit=21,file="altitude.txt",action="write")
+    DO k=1,kte!-1
+      p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top
+      altitude(k) = interp_0_log( zk, pd_in, p_level, nl_in )
+      write (21,*) altitude(k)
+    ENDDO
+    close(unit=21)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -941 +968 @@
 !---------------------------------------------------------------------------
 
-      subroutine get_sounding( zk, p, p_dry, theta, rho, &
+      subroutine get_sounding( zk, p, p_dry, theta, Tz, rho, &
                                u, v, qv, dry, nl_max, nl_in )
       implicit none
@@ -947 +974 @@
       integer nl_max, nl_in
       real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), &
-           u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max)
+           u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max),Tz(nl_max)
       logical dry
 
@@ -1112 +1139 @@
           v(k) = v_input(k)
           qv(k) = qv_input(k)
+          Tz(k) = t_therm(k)
 
         enddo
@@ -1198 +1226 @@
       implicit none
       integer n
-      real pdustq(n),pdustn(n)
+      real pdustq(n+1),pdustn(n+1)
       logical end_of_file
 
@@ -1223 +1251 @@
       end subroutine read_dust
 !!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-      subroutine read_hr(hr_sw,hr_lw,n)
+      subroutine read_hr(hr_sw,hr_lw,hr_dyn,n)
       implicit none
       integer n
-      real hr_sw(n),hr_lw(n)
+      real hr_sw(n+1),hr_lw(n+1),hr_dyn(n+1)
       logical end_of_file
 
@@ -1239 +1267 @@
       do while (.not. end_of_file)
 
-        read(11,*,end=103) hr_sw(k+1),hr_lw(k+1)
-        write(*,*) k,hr_sw(k+1),hr_lw(k+1)
+        read(11,*,end=103) hr_sw(k+1),hr_lw(k+1),hr_dyn(k+1)
+        write(*,*) k,hr_sw(k+1),hr_lw(k+1),hr_dyn(k+1)
         k = k+1
-        go to 113
+        go to 114
  103    end_of_file = .true.
  114    continue