Changeset 1398 for trunk/MESOSCALE/LMD_MM_MARS/SRC

Timestamp:

Mar 18, 2015, 11:38:45 AM (10 years ago)

Author:

aslmd

Message:

changed Venus LES ini: possibility to use custom levels, correction hypso, and consistent direct-from-file mode

File:

• trunk/MESOSCALE/LMD_MM_MARS/SRC/LES/modif_mars/module_initialize_les.F.venus (modified) (15 diffs)

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

@@ -118 +118 @@
 !!MARS
 
+      REAL :: pfu, pfd, phm
+      INTEGER :: hypsometric_opt = 1 ! classic
+      !INTEGER :: hypsometric_opt = 2 ! Wee et al. 2012 correction
+
+
+
 #ifdef DM_PARALLEL
 #    include <data_calls.inc>
@@ -164 +170 @@
 
 
-!  stretch_grid = .true.
+   stretch_grid = .true.
 !  FOR LES, set stretch to false
    stretch_grid = .false.
@@ -170 +176 @@
 !   z_scale = .50
    z_scale = .40
+   z_scale = .25
+   z_scale = .15
    pi = 2.*asin(1.0)
    write(6,*) ' pi is ',pi
@@ -232 +240 @@
    grid%step_number = 0
 
-! set up the grid
-
+!! set up the grid
+!
    IF (stretch_grid) THEN ! exponential stretch for eta (nearly constant dz)
      DO k=1, kde
@@ -240 +248 @@
      ENDDO
    ELSE
-
-!!!MARS
-grid%znw(1)=1.000
-grid%znw(2)=0.9995 !5m
-grid%znw(3)=0.9980 !20m
-grid%znw(4)=0.9950 !55m
-DO k=5, kde
-   grid%znw(k) = grid%znw(4) * ( 1. - float(k-4)/float(kde-4) )
-ENDDO
+!
 !!!!MARS
-!!
-!     DO k=1, kde
-!      grid%znw(k) = 1. - float(k-1)/float(kde-1)
-!     ENDDO
-
-   ENDIF
+!grid%znw(1)=1.000
+!grid%znw(2)=0.9995 !5m
+!grid%znw(3)=0.9980 !20m
+!grid%znw(4)=0.9950 !55m
+!DO k=5, kde
+!   grid%znw(k) = grid%znw(4) * ( 1. - float(k-4)/float(kde-4) )
+!ENDDO
+!!!!!MARS
+!!!
+!!     DO k=1, kde
+!!      grid%znw(k) = 1. - float(k-1)/float(kde-1)
+!!     ENDDO
+!
+!   ENDIF
+
+!!MARS
+!!MARS
+  open(unit=12,file='levels',form='formatted',status='old')
+  rewind(12)
+  DO k=1, kde
+  read(12,*) grid%znw(k)
+  write(6,*) 'read level ', k,grid%znw(k)
+  ENDDO
+  close(12)
+!!MARS
+!!MARS
+     ENDIF
 
    DO k=1, kde-1
@@ -419 +440 @@
 !  respect to the model eqns.
 
+   IF (hypsometric_opt == 1) THEN
     DO k  = 2,kte
       grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)*grid%mub(i,j)*grid%alb(i,k-1,j)
     ENDDO
+   ELSE IF (hypsometric_opt == 2) THEN
+    DO k = 2,kte
+      pfu = grid%mub(i,j)*grid%znw(k)   + grid%p_top
+      pfd = grid%mub(i,j)*grid%znw(k-1)   + grid%p_top
+      phm = grid%mub(i,j)*grid%znu(k-1)   + grid%p_top
+      grid%phb(i,k,j) = grid%phb(i,k-1,j) + grid%alb(i,k-1,j)*phm*LOG(pfd/pfu)
+    END DO
+   END IF
+
 
   ENDDO
@@ -501 +532 @@
 
     grid%ph_1(i,1,j) = 0.
+   IF (hypsometric_opt == 1) THEN
     DO k  = 2,kte
       grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*(       &
@@ -509 +541 @@
       grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
     ENDDO
+   ELSE IF (hypsometric_opt == 2) THEN
+
+             ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
+             ! Note that al*p approximates Rd*T and dLOG(p) does z.
+             ! Here T varies mostly linear with z, the first-order integration produces better result.
+
+               grid%ph_2(i,1,j) = grid%phb(i,1,j)
+               DO k = 2,kte
+                  pfu = grid%mu0(i,j)*grid%znw(k)   + grid%p_top
+                  pfd = grid%mu0(i,j)*grid%znw(k-1) + grid%p_top
+                  phm = grid%mu0(i,j)*grid%znu(k-1) + grid%p_top
+                  grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
+               END DO
+
+               DO k = 1,kte
+                  grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
+                  grid%ph_1(i,k,j) = grid%ph_2(i,k,j)
+               END DO
+
+   END IF
+
+
 
     IF ( wrf_dm_on_monitor() ) THEN
@@ -566 +620 @@
 !  rebalance hydrostatically
 
+   IF (hypsometric_opt == 1) THEN
+
       DO k  = 2,kte
         grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*(       &
@@ -575 +631 @@
       ENDDO
 
+   ELSE IF (hypsometric_opt == 2) THEN
+
+             ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
+             ! Note that al*p approximates Rd*T and dLOG(p) does z.
+             ! Here T varies mostly linear with z, the first-order integration produces better result.
+
+               grid%ph_2(i,1,j) = grid%phb(i,1,j)
+               DO k = 2,kte
+                  pfu = grid%mu0(i,j)*grid%znw(k)   + grid%p_top
+                  pfd = grid%mu0(i,j)*grid%znw(k-1) + grid%p_top
+                  phm = grid%mu0(i,j)*grid%znu(k-1) + grid%p_top
+                  grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
+               END DO
+
+               DO k = 1,kte
+                  grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
+                  grid%ph_1(i,k,j) = grid%ph_2(i,k,j)
+               END DO
+
+   END IF
+
+
     ENDDO
   ENDDO
@@ -609 +687 @@
       z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g
     END IF
+
     p_surf = interp_0( p_in, zk, z_at_v, nl_in )
 
@@ -879 +958 @@
       real qvf, qvf1, dz
 
+      LOGICAL :: direct_from_file = .true.
+
 !  first, read the sounding
 
       call read_sounding( p_surf, th_surf, qv_surf, &
-                          h_input, th_input, qv_input, u_input, v_input,n, nl, debug )
+                          h_input, th_input, qv_input, &
+                          u_input, v_input,n, nl, debug )
 
 ! and the therm :
 
       call read_therm(r_therm,cp_therm,p_therm,rho_therm,t_therm,n)
+
+      if(debug) write(6,*) ' number of input levels = ',nl
+      nl_in = nl
+      if(nl_in .gt. nl_max ) then
+        write(6,*) ' too many levels for input arrays ',nl_in,nl_max
+        call wrf_error_fatal ( ' too many levels for input arrays ' )
+      end if
+
+      IF (.NOT. direct_from_file) THEN
 
 ! To use r/cp as defined above, one has to recompute teta from T (default MCD computes
@@ -902 +993 @@
       ! enddo
       !endif
-
-      if(debug) write(6,*) ' number of input levels = ',nl
-
-        nl_in = nl
-        if(nl_in .gt. nl_max ) then
-          write(6,*) ' too many levels for input arrays ',nl_in,nl_max
-          call wrf_error_fatal ( ' too many levels for input arrays ' )
-        end if
 
 !  compute diagnostics,
@@ -964 +1047 @@
           qvf1 = 1.
 !!MARS 
+
             do it=1,10
               pm_input(k) = pm_input(k-1) &
@@ -982 +1066 @@
             p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g
           enddo
+
+      ELSE !IF (.NOT. direct_from_file) THEN
+        
+        do k=1,nl
+         !!!! direct input from file
+         write(6,*) '*** DIRECT INPUT FROM FILE ***'
+         pm_input(k) = p_therm(k)
+         p_input(k) = p_therm(k)
+         rho_input(k) = rho_therm(k)
+        enddo
+
+      ENDIF