Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1387 r1398  
    118118!!MARS
    119119
     120      REAL :: pfu, pfd, phm
     121      INTEGER :: hypsometric_opt = 1 ! classic
     122      !INTEGER :: hypsometric_opt = 2 ! Wee et al. 2012 correction
     123
     124
     125
    120126#ifdef DM_PARALLEL
    121127#    include <data_calls.inc>
    … …  
    164170
    165171
    166 !  stretch_grid = .true.
     172   stretch_grid = .true.
    167173!  FOR LES, set stretch to false
    168174   stretch_grid = .false.
    … …  
    170176!   z_scale = .50
    171177   z_scale = .40
     178   z_scale = .25
     179   z_scale = .15
    172180   pi = 2.*asin(1.0)
    173181   write(6,*) ' pi is ',pi
    … …  
    232240   grid%step_number = 0
    233241
    234 ! set up the grid
    235 
     242!! set up the grid
     243!
    236244   IF (stretch_grid) THEN ! exponential stretch for eta (nearly constant dz)
    237245     DO k=1, kde
    … …  
    240248     ENDDO
    241249   ELSE
    242 
    243 !!!MARS
    244 grid%znw(1)=1.000
    245 grid%znw(2)=0.9995 !5m
    246 grid%znw(3)=0.9980 !20m
    247 grid%znw(4)=0.9950 !55m
    248 DO k=5, kde
    249    grid%znw(k) = grid%znw(4) * ( 1. - float(k-4)/float(kde-4) )
    250 ENDDO
     250!
    251251!!!!MARS
    252 !!
    253 !     DO k=1, kde
    254 !      grid%znw(k) = 1. - float(k-1)/float(kde-1)
    255 !     ENDDO
    256 
    257    ENDIF
     252!grid%znw(1)=1.000
     253!grid%znw(2)=0.9995 !5m
     254!grid%znw(3)=0.9980 !20m
     255!grid%znw(4)=0.9950 !55m
     256!DO k=5, kde
     257!   grid%znw(k) = grid%znw(4) * ( 1. - float(k-4)/float(kde-4) )
     258!ENDDO
     259!!!!!MARS
     260!!!
     261!!     DO k=1, kde
     262!!      grid%znw(k) = 1. - float(k-1)/float(kde-1)
     263!!     ENDDO
     264!
     265!   ENDIF
     266
     267!!MARS
     268!!MARS
     269  open(unit=12,file='levels',form='formatted',status='old')
     270  rewind(12)
     271  DO k=1, kde
     272  read(12,*) grid%znw(k)
     273  write(6,*) 'read level ', k,grid%znw(k)
     274  ENDDO
     275  close(12)
     276!!MARS
     277!!MARS
     278     ENDIF
    258279
    259280   DO k=1, kde-1
    … …  
    419440!  respect to the model eqns.
    420441
     442   IF (hypsometric_opt == 1) THEN
    421443    DO k  = 2,kte
    422444      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)
    423445    ENDDO
     446   ELSE IF (hypsometric_opt == 2) THEN
     447    DO k = 2,kte
     448      pfu = grid%mub(i,j)*grid%znw(k)   + grid%p_top
     449      pfd = grid%mub(i,j)*grid%znw(k-1)   + grid%p_top
     450      phm = grid%mub(i,j)*grid%znu(k-1)   + grid%p_top
     451      grid%phb(i,k,j) = grid%phb(i,k-1,j) + grid%alb(i,k-1,j)*phm*LOG(pfd/pfu)
     452    END DO
     453   END IF
     454
    424455
    425456  ENDDO
    … …  
    501532
    502533    grid%ph_1(i,1,j) = 0.
     534   IF (hypsometric_opt == 1) THEN
    503535    DO k  = 2,kte
    504536      grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*(       &
    … …  
    509541      grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
    510542    ENDDO
     543   ELSE IF (hypsometric_opt == 2) THEN
     544
     545             ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
     546             ! Note that al*p approximates Rd*T and dLOG(p) does z.
     547             ! Here T varies mostly linear with z, the first-order integration produces better result.
     548
     549               grid%ph_2(i,1,j) = grid%phb(i,1,j)
     550               DO k = 2,kte
     551                  pfu = grid%mu0(i,j)*grid%znw(k)   + grid%p_top
     552                  pfd = grid%mu0(i,j)*grid%znw(k-1) + grid%p_top
     553                  phm = grid%mu0(i,j)*grid%znu(k-1) + grid%p_top
     554                  grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
     555               END DO
     556
     557               DO k = 1,kte
     558                  grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
     559                  grid%ph_1(i,k,j) = grid%ph_2(i,k,j)
     560               END DO
     561
     562   END IF
     563
     564
    511565
    512566    IF ( wrf_dm_on_monitor() ) THEN
    … …  
    566620!  rebalance hydrostatically
    567621
     622   IF (hypsometric_opt == 1) THEN
     623
    568624      DO k  = 2,kte
    569625        grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*(       &
    … …  
    575631      ENDDO
    576632
     633   ELSE IF (hypsometric_opt == 2) THEN
     634
     635             ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
     636             ! Note that al*p approximates Rd*T and dLOG(p) does z.
     637             ! Here T varies mostly linear with z, the first-order integration produces better result.
     638
     639               grid%ph_2(i,1,j) = grid%phb(i,1,j)
     640               DO k = 2,kte
     641                  pfu = grid%mu0(i,j)*grid%znw(k)   + grid%p_top
     642                  pfd = grid%mu0(i,j)*grid%znw(k-1) + grid%p_top
     643                  phm = grid%mu0(i,j)*grid%znu(k-1) + grid%p_top
     644                  grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
     645               END DO
     646
     647               DO k = 1,kte
     648                  grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
     649                  grid%ph_1(i,k,j) = grid%ph_2(i,k,j)
     650               END DO
     651
     652   END IF
     653
     654
    577655    ENDDO
    578656  ENDDO
    … …  
    609687      z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g
    610688    END IF
     689
    611690    p_surf = interp_0( p_in, zk, z_at_v, nl_in )
    612691
    … …  
    879958      real qvf, qvf1, dz
    880959
     960      LOGICAL :: direct_from_file = .true.
     961
    881962!  first, read the sounding
    882963
    883964      call read_sounding( p_surf, th_surf, qv_surf, &
    884                           h_input, th_input, qv_input, u_input, v_input,n, nl, debug )
     965                          h_input, th_input, qv_input, &
     966                          u_input, v_input,n, nl, debug )
    885967
    886968! and the therm :
    887969
    888970      call read_therm(r_therm,cp_therm,p_therm,rho_therm,t_therm,n)
     971
     972      if(debug) write(6,*) ' number of input levels = ',nl
     973      nl_in = nl
     974      if(nl_in .gt. nl_max ) then
     975        write(6,*) ' too many levels for input arrays ',nl_in,nl_max
     976        call wrf_error_fatal ( ' too many levels for input arrays ' )
     977      end if
     978
     979      IF (.NOT. direct_from_file) THEN
    889980
    890981! To use r/cp as defined above, one has to recompute teta from T (default MCD computes
    … …  
    902993      ! enddo
    903994      !endif
    904 
    905       if(debug) write(6,*) ' number of input levels = ',nl
    906 
    907         nl_in = nl
    908         if(nl_in .gt. nl_max ) then
    909           write(6,*) ' too many levels for input arrays ',nl_in,nl_max
    910           call wrf_error_fatal ( ' too many levels for input arrays ' )
    911         end if
    912995
    913996!  compute diagnostics,
    … …  
    9641047          qvf1 = 1.
    9651048!!MARS
     1049
    9661050            do it=1,10
    9671051              pm_input(k) = pm_input(k-1) &
    … …  
    9821066            p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g
    9831067          enddo
     1068
     1069      ELSE !IF (.NOT. direct_from_file) THEN
     1070       
     1071        do k=1,nl
     1072         !!!! direct input from file
     1073         write(6,*) '*** DIRECT INPUT FROM FILE ***'
     1074         pm_input(k) = p_therm(k)
     1075         p_input(k) = p_therm(k)
     1076         rho_input(k) = rho_therm(k)
     1077        enddo
     1078
     1079      ENDIF
    9841080
    9851081
Note: See TracChangeset for help on using the changeset viewer.