Ignore:
Timestamp:
Jul 21, 2017, 4:07:38 PM (7 years ago)
Author:
mlefevre
Message:

Add of Venus mesoscale structure + update of Venus LES

Location:
trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2
Files:
12 added
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM

    r1611 r1724  
    112112state  real  WSTAR     ij   misc  1  -  rd   "WSTAR"     "FREE CONVECTION VELOCITY FROM TH" "m/s"    #SAVEMARS2 wstar
    113113state  real  Z0SET     ij   misc  1  -  rd   "Z0SET"     "SET SURFACE ROUGHNESS"           "m"       #SAVEMARS2 z0
     114state  real  HR_DYN     ikj   misc  1  -  -  "HR_DYN"     "HEATING RATE DYN"                 "K/s"
     115state  real  DDT     ikj   misc  1  -  -  "DDT"     "TEMPERATURE TENDENCY"                 "K/s"
     116state  real  DT_RAD     ikj   misc  1  -  -  "DT_RAD"     "RADIATIVE TEMPERATURE TENDENCY"                 "K/s"
     117state  real  DT_VDF     ikj   misc  1  -  -  "DT_VDF"     "PBL CLMAIN TEMPERATURE TENDENCY"                 "K/s"
     118state  real  DT_AJS     ikj   misc  1  -  -  "DT_AJS"     "AJUSTEMENT SEC TEMPERATURE TENDENCY"                 "K/s"
    114119
    115120####

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM.newphys

    r1635 r1724  
    8080state  real  HR_SW     ikj   misc  1  -  h  "HR_SW"     "HEATING RATE SW"                 "K/s"
    8181state  real  HR_LW     ikj   misc  1  -  h  "HR_LW"     "HEATING RATE LW"                 "K/s"
     82state  real  HR_DYN     ikj   misc  1  -  -  "HR_DYN"     "HEATING RATE DYN"                 "K/s"
     83state  real  DDT     ikj   misc  1  -  -  "DDT"     "TEMPERATURE TENDENCY"                 "K/s"
     84state  real  DT_RAD     ikj   misc  1  -  -  "DT_RAD"     "RADIATIVE TEMPERATURE TENDENCY"                 "K/s"
     85state  real  DT_VDF     ikj   misc  1  -  -  "DT_VDF"     "PBL CLMAIN TEMPERATURE TENDENCY"                 "K/s"
     86state  real  DT_AJS     ikj   misc  1  -  -  "DT_AJS"     "AJUSTEMENT SEC TEMPERATURE TENDENCY"                 "K/s"
    8287state  real  SWDOWNZ    ij   misc  1  -  h  "SWDOWNZ"   "DOWNWARD SW FLUX AT SURFACE"     "W m-2"
    8388#### DUST CYCLE
    … …  
    106111#### - effacer le h si l'on veut eviter les sorties
    107112####
    108 state  real  -          ikjftb  scalar  1  -  -                         -
    109 state  real  qh2o       ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QH2O"          "Water vapor mixing ratio"    "kg kg-1"
    110 state  real  qh2o_ice   ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QH2O_ICE"      "Water ice mixing ratio"      "kg kg-1"
    111 state  real  qdust      ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QDUST"         "Dust mixing ratio"           "kg kg-1"
    112 state  real  qdustn     ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QDUSTN"        "Dust_number mixing ratio"    "kg kg-1"
    113 state  real  qccn       ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QCCN"          "CCN mixing ratio"            "kg kg-1"
    114 state  real  qccnn      ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QCCNN"         "CCN number mixing ratio"     "kg kg-1"
    115 state  real  qco2       ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QCO2"          "CO2 mixing ratio"            "kg kg-1"
    116 state  real  qco2_ice   ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QCO2_ICE"      "CO2 ice mixing ratio"        "kg kg-1"
    117 state  real  qccn_co2   ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QCCN_CO2"      "CO2 CCN mixing ratio"        "kg kg-1"
    118 state  real  qccnn_co2  ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QCCNN_CO2"     "CO2 CCN number mixing ratio" "kg kg-1"
     113state  real  -         ikjftb  scalar  1  -  -                         -
     114state  real  qh2o      ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QH2O"          "Water vapor mixing ratio"   "kg kg-1"
     115state  real  qh2o_ice  ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QH2O_ICE"      "Water ice mixing ratio"     "kg kg-1"
     116state  real  qdust     ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QDUST"         "Dust mixing ratio"          "kg kg-1"
     117state  real  qdustn    ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QDUSTN"        "Dust_number mixing ratio"   "kg kg-1"
     118state  real  qccn      ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QCCN"         "CCN mixing ratio"          "kg kg-1"
     119state  real  qccnn     ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QCCNN"        "CCN number mixing ratio"   "kg kg-1"
     120state  real  qco2      ikjftb  scalar  1  -  i01rhusdf=(bdy_interp:dt)  "QCO2"          "CO2 mixing ratio"           "kg kg-1"
    119121state  real  chem_co   ikjftb  scalar  1  -  i01rusdf=(bdy_interp:dt)  "chem_co"     ""   ""
    120122state  real  chem_o    ikjftb  scalar  1  -  i01rusdf=(bdy_interp:dt)  "chem_o"      ""   ""
    … …  
    130132state  real  chem_n2   ikjftb  scalar  1  -  i01rusdf=(bdy_interp:dt)  "chem_n2"     ""   ""
    131133state  real  chem_ar   ikjftb  scalar  1  -  i01rusdf=(bdy_interp:dt)  "chem_ar"     ""   ""
     134####
    132135####
    133136####MARS MARS MARS MARS MARS
    … …  
    12671270rconfig   logical pert_coriolis           namelist,dynamics     max_domains  .false.  irh  "pert_coriolis"    ""   ""
    12681271rconfig   logical mix_full_fields         namelist,dynamics     max_domains  .false.  irh  "mix_full_field"   ""   ""
    1269 rconfig   real    base_pres               namelist,dynamics     1             610.     h    "base_pres"  "Mars Base state pressure, real only"      "Pa"
    1270 rconfig   real    base_temp               namelist,dynamics     1             220.     h    "base_temp"  "Mars Base state sea level temperature, real only"      "K"
    1271 rconfig   real    base_lapse              namelist,dynamics     1              30.     h    "base_lapse" "Mars Base state temperature difference between base pres and 1/e of atm depth - do not change, real only"      "K"
     1272rconfig   real    base_pres               namelist,dynamics     1             9.2e6     h    "base_pres"  "Mars Base state pressure, real only"      "Pa"
     1273rconfig   real    base_temp               namelist,dynamics     1             735.     h    "base_temp"  "Mars Base state sea level temperature, real only"      "K"
     1274rconfig   real    base_lapse              namelist,dynamics     1             100.     h    "base_lapse" "Mars Base state temperature difference between base pres and 1/e of atm depth - do not change, real only"      "K"
    12721275rconfig   real    tiso                    namelist,dynamics     1             140.     h    "tiso"    "Mars isothermal value for high alt, real only"  "K"
    12731276
    … …  
    13541357package   newwater     mars==11                     -              scalar:qh2o,qh2o_ice,qdust,qdustn
    13551358package   newwatersca  mars==12                     -              scalar:qh2o,qh2o_ice,qdust,qdustn,qccn,qccnn
    1356 package   cloudco2     mars==32                     -              scalar:qh2o,qh2o_ice,qdust,qdustn,qccn,qccnn,qco2,qco2_ice,qccn_co2,qccnn_co2
    13571359package   photochem    mars==42                     -              scalar:qco2,chem_co,chem_o,chem_o1d,chem_o2,chem_o3,chem_h,chem_h2,chem_oh,chem_ho2,chem_h2o2,chem_ch4,chem_n2,chem_ar,qh2o_ice,qh2o,qdust,qdustn
    13581360##### MARS OPTIONS

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_big_step_utilities_em.F

    r11 r1724  
    1818   USE module_configure
    1919   USE module_wrf_error
    20 
    2120CONTAINS
    2221
    … …  
    569568      ENDDO
    570569      ENDDO
    571 
     570      !print*,'muu(i,j)',muu(1,:)
    572571   ELSE IF (name .EQ. 'v')THEN
    573572
    … …  
    588587      ENDDO
    589588      ENDDO
    590 
     589      !print*,'muv(i,j)',muv(1,:)
    591590   ELSE IF (name .EQ. 'w')THEN
    592591      itf=MIN(ite,ide-1)
    … …  
    621620      ENDDO
    622621      ENDDO
    623    
     622      !print*,'field(i,k,j)',field(1,:,1)
     623      !print*,'(mu(i,j)+mub(i,j))',(mu(1,:)+mub(1,:))   
    624624   ENDIF
    625625
    … …  
    10521052
    10531053  IF (non_hydrostatic) THEN
    1054 
     1054      print*,'n_moist',n_moist
     1055      print*,'PARAM_FIRST_SCALAR',PARAM_FIRST_SCALAR
    10551056      IF (n_moist >= PARAM_FIRST_SCALAR ) THEN 
    1056 
     1057        print*,'p(i,k,j)_before',p(1,:,1)
    10571058        DO j=jts,jtf
    10581059        DO k=kts,ktf
    … …  
    10751076        ENDDO
    10761077        ENDDO
    1077 
     1078        print*,'p(i,k,j)_after',p(1,:,1)
    10781079      ELSE
    10791080
    … …  
    40754076   INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
    40764077   INTEGER :: i, j, k
    4077    REAL    :: w1, w2, z0, z1, z2
     4078   REAL    :: w1, w2, z0, z1, z2,temp
    40784079
    40794080!-----------------------------------------------------------------------
    … …  
    40994100
    41004101!  compute thermodynamics and velocities at pressure points
    4101 
    41024102    do j = j_start,j_end
    41034103    do k = k_start, k_end
    … …  
    41074107      p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
    41084108      pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
    4109       t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
     4109      IF (planet .eq. "venus" ) THEN
     4110        t_phy(i,k,j)= (th_phy(i,k,j)**nu - nu*(TT00**nu)*log((p1000mb/p_phy(i,k,j))**rcp))**(1/nu)
     4111      ELSE
     4112        t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
     4113      ENDIF
    41104114      rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
    41114115      mu_3d(i,k,j) = mu(i,j)

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_init_utilities.F

    r11 r1724  
    9191 END FUNCTION interp_0
    9292
     93 real function interp_0_log( v_in,  &
     94                         p_in, p_out, nz_in  )
     95 implicit none
     96 integer nz_in, nz_out
     97 real    v_in(nz_in), p_in(nz_in)
     98 real    p_out
     99
     100 integer kp, k, im, ip
     101 logical interp, increasing_z
     102 real    height, w1, w2
     103 logical debug
     104 parameter ( debug = .false. )
     105
     106! does vertical coordinate increase or decrease with increasing k?
     107! set offset appropriately
     108
     109 height = p_out
     110
     111 if(debug) write(6,*) ' height in interp_0 ',height
     112
     113 if (p_in(nz_in) .gt. p_in(1)) then
     114
     115    if(debug) write(6,*) ' monotonic increase in z in interp_0 '
     116    IF (height > p_in(nz_in)) then
     117      if(debug) write(6,*) ' point 1 in interp_0 '
     118      w2 = log(p_in(nz_in)/height)/log(p_in(nz_in)/p_in(nz_in-1))
     119      w1 = 1.-w2
     120      interp_0_log = w1*v_in(nz_in) + w2*v_in(nz_in-1)
     121    ELSE IF (height < p_in(1)) then
     122      if(debug) write(6,*) ' point 2 in interp_0 '
     123      w2 = log(p_in(2)/height)/log(p_in(2)/p_in(1))
     124      w1 = 1.-w2
     125      interp_0_log = w1*v_in(2) + w2*v_in(1)
     126    ELSE
     127      if(debug) write(6,*) ' point 3 in interp_0 '
     128      interp = .false.
     129      kp = nz_in
     130      DO WHILE ( (interp .eqv. .false.) .and. (kp .ge. 2) )
     131        IF(   ((p_in(kp)   .ge. height) .and.     &
     132               (p_in(kp-1) .le. height))        )   THEN
     133          w2 = log(height/p_in(kp))/log(p_in(kp-1)/p_in(kp))
     134          w1 = 1.-w2
     135          interp_0_log = w1*v_in(kp) + w2*v_in(kp-1)
     136          if(debug) write(6,*) ' interp data, kp, w1, w2 ',kp, w1, w2
     137          if(debug) write(6,*) ' interp data, v_in(kp), v_in(kp-1), interp_0_p ', &
     138                     v_in(kp), v_in(kp-1), interp_0_log
     139          interp = .true.
     140        END IF
     141        kp = kp-1
     142      ENDDO
     143    ENDIF
     144
     145 else
     146
     147    if(debug) write(6,*) ' monotonic decrease in z in interp_0 '
     148
     149    IF (height < p_in(nz_in)) then
     150      if(debug) write(6,*) ' point 1 in interp_0 '
     151      w2 = log(p_in(nz_in)/height)/log(p_in(nz_in)/p_in(nz_in-1))
     152      w1 = 1.-w2
     153      interp_0_log = w1*v_in(nz_in) + w2*v_in(nz_in-1)
     154    ELSE IF (height > p_in(1)) then
     155      if(debug) write(6,*) ' point 2 in interp_0 '
     156      w2 = log(p_in(2)/height)/log(p_in(2)/p_in(1))
     157      w1 = 1.-w2
     158      interp_0_log = w1*v_in(2) + w2*v_in(1)
     159    ELSE
     160      if(debug) write(6,*) ' point 3 in interp_0 '
     161      interp = .false.
     162      kp = nz_in
     163      height = p_out
     164      DO WHILE ( (interp .eqv. .false.) .and. (kp .ge. 2) )
     165        IF(   ((p_in(kp)   .le. height) .and.     &
     166               (p_in(kp-1) .ge. height))             )   THEN
     167          w2 = log(height/p_in(kp))/log(p_in(kp-1)/p_in(kp))
     168          w1 = 1.-w2
     169          interp_0_log = w1*v_in(kp) + w2*v_in(kp-1)
     170          interp = .true.
     171        END IF
     172        kp = kp-1
     173      ENDDO
     174    ENDIF
     175
     176 end if
     177
     178 return
     179 END FUNCTION interp_0_log
     180
    93181END MODULE module_init_utilities
    94182

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_real.F

    r1662 r1724  
    2121   USE module_soil_pre
    2222   USE module_date_time
     23   use module_init_utilities
    2324#ifdef DM_PARALLEL
    2425   USE module_dm
    … …  
    118119      REAL :: max_dz
    119120
    120 !      INTEGER , PARAMETER :: nl_max = 1000
     121      INTEGER , PARAMETER :: nl_max = 1000
     122      INTEGER :: nl_in
    121123!      REAL , DIMENSION(nl_max) :: grid%em_dn
    122 
     124      REAL, DIMENSION(nl_max) :: pd_in
     125      REAL, DIMENSION(nl_max) :: profdustq,profdustn
    123126integer::oops1,oops2
    124127
    … …  
    149152      !LOGICAL :: interp_theta = .false. ! Wee et al. 2012 correction
    150153      REAL :: pfu, pfd, phm
    151 
     154      REAL :: tpot
    152155
    153156#ifdef DM_PARALLEL
    … …  
    10801083!!END DO
    10811084
    1082 
     1085IF ( planet == "mars" ) then
    10831086!--get vertical size of the GCM input array and allocate new stuff
    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 )
     1087  sizegcm=SIZE(grid%em_rh_gc(1,:,1))
     1088  ALLOCATE(sig(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1))
     1089  !ALLOCATE(ap(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1))
     1090  ALLOCATE(bp(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1))
     1091
     1092  DO j = jts , MIN ( jde-1 , jte )
     1093  DO i = its , MIN (ide-1 , ite )
    10911094
    10921095!!! Define old sigma levels for each column
    1093         sig(i,:,j)=grid%em_p_gc(i,:,j)/grid%em_psfc_gc(i,j)
     1096          sig(i,:,j)=grid%em_p_gc(i,:,j)/grid%em_psfc_gc(i,j)
    10941097
    10951098!!! Compute new sigma levels from old sigma levels with GCM (low-res) and WRF (hi-res) surface pressure
    10961099!!!                        (dimlevs,sigma_gcm, ps_gcm,                ps_hr,         sigma_hr)
    1097         CALL build_sigma_hr(sizegcm,sig(i,:,j),grid%em_psfc_gc(i,j),grid%psfc(i,j),bp(i,:,j))
     1100          CALL build_sigma_hr(sizegcm,sig(i,:,j),grid%em_psfc_gc(i,j),grid%psfc(i,j),bp(i,:,j))
    10981101
    10991102!!! Calculate new pressure levels
    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)
     1103          grid%em_pd_gc(i,:,j)=bp(i,:,j)*grid%psfc(i,j)
     1104
     1105  END DO
     1106  END DO
     1107
     1108  DEALLOCATE(sig)
     1109  DEALLOCATE(bp)
    11071110
    11081111!!****MARS who knows...
    1109 grid%em_rh_gc(:,:,:)=0.
     1112  grid%em_rh_gc(:,:,:)=0.
    11101113
    11111114
    … …  
    11131116!grid%em_pd_gc=grid%em_p_gc
    11141117!!****MARS
    1115 
     1118ELSE ! VENUS
    11161119
    11171120
    … …  
    11201123         !!  dry top pressure (constant).
    11211124         !
    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 
     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 )
     1129ENDIF
     1130IF ( planet == "mars" ) then
    11271131!!****MARS
    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
     1132  DO j = jts , MIN ( jde-1 , jte )
     1133  DO i = its , MIN (ide-1 , ite )
     1134
     1135     grid%em_mu0(i,j) = grid%psfc(i,j) - grid%p_top
     1136
     1137  END DO
     1138  END DO
    11351139!!****MARS
    1136 
     1140ELSE ! VENUS
    11371141   
    11381142         !!  Compute the dry, hydrostatic surface pressure.
    11391143         !
    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 )
     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 )
     1148ENDIF
    11441149!!****MARS: voir remarques dans la routine
    11451150!!****MARS: dry hydrostatic pressure comes from the GCM ...
    … …  
    22942299IF (( i .EQ. its ) .AND. ( j .EQ. jts )) print *, temp, k
    22952300!!! MODIF WRFV3.1 - parameter tiso
    2296                grid%em_t_init(i,k,j) = temp*(p00/grid%em_pb(i,k,j))**(r_d/cp) - t0
     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
    22972306               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
    22982307            END DO
    … …  
    26012610!!NB: q2 is used for other purpose ...
    26022611
    2603 
    26042612!      END IF
    26052613
    … …  
    27802788temp1 = MAX(tiso,t00+A*LOG(grid%em_pb(i,k,j)/p00))
    27812789temp2 = MAX(tiso,t00+A*LOG(           pb_int/p00))
    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
     2790IF (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
     2793ELSE
     2794ENDIF
    27842795               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
    27852796            END DO
    … …  
    49694980!            temp =             t00 + A*LOG(pb/p00)
    49704981temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
    4971             t_init = temp*(p00/pb)**(r_d/cp) - t0
     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
    49724987            alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
    49734988         END DO
    … …  
    50055020!            temp =             t00 + A*LOG(pb/p00)
    50065021temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
    5007             t_init = temp*(p00/pb)**(r_d/cp) - t0
     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
    50085027            alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
    50095028            znw(k+1) = znw(k) - dz*g / ( mub*alb(k) )
    … …  
    50255044!                  temp =             t00 + A*LOG(pb/p00)
    50265045temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
    5027                   t_init = temp*(p00/pb)**(r_d/cp) - t0
     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
    50285051                  alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
    50295052                  znw(k+1) = znw(k) - dz*g / ( mub*alb(k) )
    … …  
    50525075!               temp =             t00 + A*LOG(pb/p00)
    50535076temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
    5054                t_init = temp*(p00/pb)**(r_d/cp) - t0
     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
    50555082               alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
    50565083            END DO
    … …  
    57205747      end subroutine build_sigma_hr
    57215748
    5722 
    5723 
    5724 
    5725 
    57265749END MODULE module_initialize
    57275750

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/solve_em.F

    r1590 r1724  
    160160  ijds = min(ids, jds)
    161161  ijde = max(ide, jde)
    162 
     162  !num_moist = 3
    163163  num_3d_m        = num_moist
    164164  num_3d_c        = num_chem
    … …  
    892892     &       ,HR_SW=grid%HR_SW   &
    893893     &       ,HR_LW=grid%HR_LW   &
     894     &       ,HR_DYN=grid%HR_DYN   &
     895     &       ,DDT=grid%DDT  &
     896     &       ,DT_RAD=grid%DT_RAD  &
     897     &       ,DT_VDF=grid%DT_VDF  &
     898     &       ,DT_AJS=grid%DT_AJS  &
    894899     &       ,SWDOWNZ=grid%SWDOWNZ   &
    895900     &       ,TAU_DUST=grid%TAU_DUST   &

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_mars_lmd_new/update_outputs_physiq_mod.F

    r1590 r1724  
    160160            SWDOWNZ,TAU_DUST,QSURFDUST,&
    161161            MTOT,ICETOT,TAU_ICE,&
    162             HR_SW,HR_LW,&
     162            HR_SW,HR_LW,HR_DYN,DT,DTRAD,DTVDF,DTAJS,&
    163163            RDUST,VMR_ICE,RICE)
    164164
    … …  
    171171     MTOT,ICETOT,TAU_ICE
    172172   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: &
    173      HR_SW,HR_LW,RDUST,VMR_ICE,RICE
     173     HR_SW,HR_LW,RDUST,VMR_ICE,RICE,HR_DYN,DT,DTRAD,DTVDF,DTAJS
    174174   INTEGER :: i,j,subs
    175175

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_prescribed_void_lmd_new/update_outputs_physiq_mod.F

    r1603 r1724  
    5252SUBROUTINE update_outputs_physiq_turb( &
    5353            ims,ime,jms,jme,kms,kme,&
    54             ips,ipe,jps,jpe,&
     54            ips,ipe,jps,jpe,kps,kpe,&
    5555            M_Q2,M_WSTAR,&
    5656            HFMAX,ZMAX,USTM,HFX)
    … …  
    6767!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    6868SUBROUTINE update_outputs_physiq_diag( &
    69             ims,ime,jms,jme,kms,kme,&
    70             ips,ipe,jps,jpe,kps,kpe,&
    71             SWDOWNZ,TAU_DUST,QSURFDUST,&
    72             MTOT,ICETOT,TAU_ICE,&
    73             HR_SW,HR_LW,&
    74             RDUST,VMR_ICE,RICE)
     69           ims,ime,jms,jme,kms,kme,&
     70           ips,ipe,jps,jpe,kps,kpe,&
     71           SWDOWNZ,TAU_DUST,QSURFDUST,&
     72           MTOT,ICETOT,TAU_ICE,&
     73           HR_SW,HR_LW,HR_DYN,DT,DTRAD,DTVDF,DTAJS,&
     74           RDUST,VMR_ICE,RICE)
    7575
    7676   INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme
    … …  
    8080     MTOT,ICETOT,TAU_ICE
    8181   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: &
    82      HR_SW,HR_LW,RDUST,VMR_ICE,RICE
     82     HR_SW,HR_LW,HR_DYN,RDUST,VMR_ICE,RICE,DT,DTRA,DTVDF,DTAJS
    8383   INTEGER :: i,j,k,subs
    84    REAL, DIMENSION(kms:kme) :: phr_sw,phr_lw
     84   REAL, DIMENSION(kms:kme) :: phr_sw,phr_lw,phr_dyn
    8585
    8686  open(17,file='prescribed_sw.txt',form='formatted',status='old')
    8787  open(18,file='prescribed_lw.txt',form='formatted',status='old')
     88  open(19,file='prescribed_dyn.txt',form='formatted',status='old')
    8889  rewind(17)
    8990  rewind(18)
     91  rewind(19)
    9092  DO k=1,kpe-kps+1
    9193    read(17,*) phr_sw(k)
    9294    read(18,*) phr_lw(k)
     95    read(19,*) phr_dyn(k)
    9396  ENDDO
     97  close(19)
     98  close(18)
    9499  close(17)
    95   close(18)
    96   print*,'phr_lw',phr_lw
     100 !print*,'phr_lw',phr_lw
    97101  DO j=jps,jpe
    98102  DO i=ips,ipe
    99103    HR_SW(i,:,j)=phr_sw(:)
    100104    HR_LW(i,:,j)=phr_lw(:)
     105    HR_DYN(i,:,j)=phr_dyn(:)
    101106  ENDDO
    102107  ENDDO

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_venus_lmd_new/callphysiq_mod.F

    r1634 r1724  
    1212CONTAINS
    1313
    14 SUBROUTINE call_physiq(planet_type, klon,llm,nqtot,                       &
     14SUBROUTINE call_physiq(planet_type, klon,llm,nqtot,tname,                 &
    1515                       debut_split,lafin_split,                           &
    1616                       jD_cur,jH_cur_split,zdt_split,                     &
    … …  
    2323
    2424  USE physiq_mod, ONLY: physiq
    25   USE module_model_constants, only : p0,rcp
     25  USE module_model_constants, only : p0,rcp,cp
     26  use moyzon_mod, only: tmoy
    2627  IMPLICIT NONE
    2728
    … …  
    3132  INTEGER,INTENT(IN) :: llm  ! number of atmospheric layers
    3233  INTEGER,INTENT(IN) :: nqtot ! number of tracers
     34  CHARACTER(len=*),INTENT(IN) :: tname(nqtot) ! tracer names
    3335  LOGICAL,INTENT(IN) :: debut_split ! .true. if very first call to physics
    3436  LOGICAL,INTENT(IN) :: lafin_split ! .true. if last call to physics
    35   REAL*4,INTENT(IN) :: JD_cur ! Julian day
    36   REAL*4,INTENT(IN) :: JH_cur_split ! Julian hour (fraction of day)
    37   REAL*4,INTENT(IN) :: zdt_split ! time step over which the physics are evaluated
    38   REAL*4,INTENT(IN) :: zplev_omp(klon,llm+1) ! interlayer pressure (Pa)
    39   REAL*4,INTENT(IN) :: zplay_omp(klon,llm) ! mid-layer pressure (Pa)
    40   REAL*4,INTENT(INOUT) :: zpk_omp(klon,llm)
    41   REAL*4,INTENT(IN) :: zphi_omp(klon,llm) ! geopotential at midlayer
    42   REAL*4,INTENT(INOUT) :: zphis_omp(klon) ! surface geopotential
    43   REAL*4,INTENT(INOUT) :: presnivs_omp(llm) ! approximate pressure of atm. layers
    44   REAL*4,INTENT(IN) :: zufi_omp(klon,llm) ! zonal wind (m/s)
    45   REAL*4,INTENT(IN) :: zvfi_omp(klon,llm) ! meridional wind (m/s)
    46   REAL*4,INTENT(INOUT) :: zrfi_omp(klon,llm) ! relative wind vorticity, in s-1
    47   REAL*4,INTENT(IN) :: ztfi_omp(klon,llm) ! temperature (K)
    48   REAL*4,INTENT(IN) :: zqfi_omp(klon,llm,nqtot) ! tracers (*/kg of air)
    49   REAL*4,INTENT(IN) :: flxwfi_omp(klon,llm) ! Vertical mass flux on lower mesh interfaces (kg/s)
     37  REAL*8,INTENT(IN) :: JD_cur ! Julian day
     38  REAL*8,INTENT(IN) :: JH_cur_split ! Julian hour (fraction of day)
     39  REAL*8,INTENT(IN) :: zdt_split ! time step over which the physics are evaluated
     40  REAL*8,INTENT(IN) :: zplev_omp(klon,llm+1) ! interlayer pressure (Pa)
     41  REAL*8,INTENT(IN) :: zplay_omp(klon,llm) ! mid-layer pressure (Pa)
     42  REAL*8,INTENT(INOUT) :: zpk_omp(klon,llm)
     43  REAL*8,INTENT(IN) :: zphi_omp(klon,llm) ! geopotential at midlayer
     44  REAL*8,INTENT(INOUT) :: zphis_omp(klon) ! surface geopotential
     45  REAL*8,INTENT(INOUT) :: presnivs_omp(llm) ! approximate pressure of atm. layers
     46  REAL*8,INTENT(IN) :: zufi_omp(klon,llm) ! zonal wind (m/s)
     47  REAL*8,INTENT(IN) :: zvfi_omp(klon,llm) ! meridional wind (m/s)
     48  REAL*8,INTENT(INOUT) :: zrfi_omp(klon,llm) ! relative wind vorticity, in s-1
     49  REAL*8,INTENT(IN) :: ztfi_omp(klon,llm) ! temperature (K)
     50  REAL*8,INTENT(IN) :: zqfi_omp(klon,llm,nqtot) ! tracers (*/kg of air)
     51  REAL*8,INTENT(IN) :: flxwfi_omp(klon,llm) ! Vertical mass flux on lower mesh interfaces (kg/s)
    5052  ! tendencies (in */s) from the physics:
    51   REAL*4,INTENT(OUT) :: zdufi_omp(klon,llm) ! tendency on zonal winds
    52   REAL*4,INTENT(OUT) :: zdvfi_omp(klon,llm) ! tendency on meridional winds
    53   REAL*4,INTENT(OUT) :: zdtfi_omp(klon,llm) ! tendency on temperature
    54   REAL*4,INTENT(OUT) :: zdqfi_omp(klon,llm,nqtot) ! tendency on tracers
    55   REAL*4,INTENT(OUT) :: zdpsrf_omp(klon) ! tendency on surface pressure
     53  REAL*8,INTENT(OUT) :: zdufi_omp(klon,llm) ! tendency on zonal winds
     54  REAL*8,INTENT(OUT) :: zdvfi_omp(klon,llm) ! tendency on meridional winds
     55  REAL*8,INTENT(OUT) :: zdtfi_omp(klon,llm) ! tendency on temperature
     56  REAL*8,INTENT(OUT) :: zdqfi_omp(klon,llm,nqtot) ! tendency on tracers
     57  REAL*8,INTENT(OUT) :: zdpsrf_omp(klon) ! tendency on surface pressure
     58  REAL*8 :: zplevmoy(llm+1) ! planet-averaged mean pressure (Pa) at interfaces
     59  REAL*8 :: ztmoy(llm)
    5660
    5761!  ! Local variables
    … …  
    7781
    7882! Set dummy variables for Mars to zero (additional and prob useless security)
    79   zpk_omp=(zplay_omp/p0)**rcp
    80   zphis_omp(:)=zphi_omp(1,1)
     83! NB: tname already filled with tracers' names (though not used here)
     84  zpk_omp(1:klon,1:llm)=cp*((zplay_omp(1:klon,1:llm)/p0)**rcp)
     85  !print*,'zpk_omp',zpk_omp(1,:)
     86  zphis_omp(1:klon)=zphi_omp(1:klon,1)
    8187  presnivs_omp(:)=0.
    8288  zrfi_omp(:,:)=0.
     89  ztmoy(:)=ztfi_omp(1,:)
     90  zplevmoy(:)=zplev_omp(1,:)
    8391! Call physics package with required inputs/outputs
    8492  CALL physiq(klon,           &
    … …  
    101109              zqfi_omp,       &
    102110              flxwfi_omp,     &
     111              zplevmoy,       &
     112              ztmoy,          &
    103113              zdufi_omp,      &
    104114              zdvfi_omp,      &
    … …  
    106116              zdqfi_omp,      &
    107117              zdpsrf_omp)
    108 
     118   deallocate(tmoy)
    109119END SUBROUTINE call_physiq
    110120

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_venus_lmd_new/iniphysiq_mod.F

    r1596 r1724  
    1010!                            dtphys,daysec,day_ini,hour_ini
    1111use comconst_mod, only : cpp
    12 use cpdet_mod, only: cpdet,ini_cpdet
     12use cpdet_phy_mod, only: init_cpdet_phy
    1313USE control_mod, ONLY: planet_type
    1414!use surfdat_h, only: ini_surfdat_h
    … …  
    1818use dimphy, only: init_dimphy
    1919USE phys_state_var_mod
    20 
     20use module_model_constants, only : nu, TT00
    2121implicit none
    2222
    … …  
    2727REAL,intent(in) :: punjours
    2828!DOUBLE PRECISION,intent(in) :: ptimestep
    29 REAL,intent(in) :: phour_ini
     29REAL*8,intent(in) :: phour_ini
    3030
    3131!real,intent(in) :: prad ! radius of the planet (m)
    … …  
    3535!real,intent(in) :: punjours ! length (in s) of a standard day [daysec]
    3636integer,intent(in) :: pdayref ! reference day of for the simulation [day_ini]
    37 real*4,intent(in) :: ptimestep !physics time step (s) [dtphys]
     37real*8,intent(in) :: ptimestep !physics time step (s) [dtphys]
    3838integer,intent(in) :: iflag_phys ! type of physics to be called
    3939
    … …  
    4343!real,intent(in) :: phour_ini   ! start time (fraction of day) of the run 0=<phour_ini<1
    4444real,intent(in) :: piphysiq   ! call physics every piphysiq dynamical timesteps
    45 
     45real*8 :: CPPP,nuvenus,T0venus
    4646
    4747! copy some fundamental parameters to physics
    … …  
    8181      cpp=pcpp
    8282      planet_type="venus"
    83       call ini_cpdet
     83      CPPP=pcpp
     84      nuvenus=nu
     85      T0venus=TT00
     86      call init_cpdet_phy(CPPP,nuvenus,T0venus)
    8487
    8588!! a few time constants initialization

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_venus_lmd_new/update_inputs_physiq_mod.F

    r1634 r1724  
    11MODULE update_inputs_physiq_mod
    2 
    3 IMPLICIT NONE
    4 
    5 CHARACTER(len=20),save,allocatable,dimension(:) :: traceurs ! tracer names
    62
    73CONTAINS
    … …  
    1713!SUBROUTINE update_inputs_physiq_slope
    1814
     15!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     16
     17
     18!c  Returns solar longitude, Ls (in deg.), from day number (in sol),
     19!c  where sol=0=Ls=0 at the northern hemisphere spring equinox
     20
     21
     22!!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1923
    2024!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    3135  REAL, INTENT(IN) :: GMT,elaps,lon_input,ls_input,lct_input
    3236  REAL,INTENT(OUT) :: MY
    33   REAL,INTENT(OUT) :: ptime,pday
     37  REAL*8,INTENT(OUT) :: ptime,pday
    3438
    3539    !
    3640    ! specified
    3741    !
    38     ptime = (GMT + elaps/3700.) !! universal time (0<ptime<1): ptime=0.5 at 12:00 UT
     42  IF (JULYR .ne. 9999) THEN
     43    ptime = (GMT + elaps/420000.) !! universal time (0<ptime<1): ptime=0.5 at 12:00 UT
    3944    ptime = MODULO(ptime,24.)   !! the two arguments of MODULO must be of the same type
    4045    ptime = ptime / 24.
    41     pday = (JULDAY - 1 + INT((3700*GMT+elaps)/88800))
     46    pday = (JULDAY - 1 + INT((420000.0*GMT+elaps)/1.008e7))
     47    pday = MODULO(int(pday),2)
     48    MY = (JULYR-2000) + (1.008e7*(JULDAY - 1)+420000.0*GMT+elaps)/2.016e7
     49    MY = INT(MY)
     50  ELSE
     51    ptime = lct_input - lon_input / 15. + elaps/(4200.)
     52    ptime = MODULO(ptime,2808.)
     53    ptime = ptime / 2808.
     54    print*,'ptime',ptime
    4255    pday = MODULO(int(pday),669)
    43     MY = (JULYR-2000) + (88800.*(JULDAY - 1)+3700.*GMT+elaps)/59496000.
    44     MY = INT(MY)
     56    MY = 2024
     57 ENDIF
    4558
    4659END SUBROUTINE update_inputs_physiq_time
    … …  
    4861!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    4962!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    50 SUBROUTINE update_inputs_physiq_tracers(nq,MARS_MODE)
     63SUBROUTINE update_inputs_physiq_tracers(nq,MARS_MODE,tname)
    5164
    5265  !use infotrac
    … …  
    5467  implicit none
    5568  INTEGER, INTENT(IN) :: nq,MARS_MODE
    56 
    57   ALLOCATE(traceurs(nq))
    58   traceurs(:)='zolbxs'
     69  CHARACTER(len=*), INTENT(INOUT) :: tname(nq) ! tracer names
     70
     71  tname(:)='zolbxs'
    5972  !call infotrac_init ! Need traceur.def
    6073
    … …  
    93106            ips,ipe,jps,jpe,&
    94107            JULYR,ngrid,nlayer,&
    95             DX,DY,MSFT,&
     108            DDX,DDY,MSFT,&
    96109            lat_input, lon_input,&
    97110            XLAT,XLONG)
    … …  
    105118   USE geometry_mod, ONLY: latitude,latitude_deg,&
    106119                           longitude,longitude_deg,&
    107                            cell_area
     120                           cell_area,dx,dy
    108121
    109122   implicit none
    … …  
    113126   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: &
    114127     MSFT,XLAT,XLONG
    115    REAL, INTENT(IN) :: dx,dy
     128   REAL, INTENT(IN) :: ddx,ddy
    116129   REAL, INTENT(IN) :: lat_input, lon_input
    117130   INTEGER :: i,j,subs
    … …  
    130143    !----------------------------------------!
    131144    !parea(subs) = dx*dy                           !! 1. idealized cases - computational grid
    132     parea(subs) = (dx/msft(i,j))*(dy/msft(i,j))    !! 2. WRF map scale factors - assume that msfx=msfy (msf=covariance)
     145    parea(subs) = (ddx/msft(i,j))*(ddy/msft(i,j))    !! 2. WRF map scale factors - assume that msfx=msfy (msf=covariance)
    133146    !parea(subs)=dx*dy/msfu(i,j)                   !! 3. special for Mercator GCM-like simulations
    134147
    … …  
    167180   IF(.not.ALLOCATED(latitude_deg)) ALLOCATE(latitude_deg(ngrid))
    168181   IF(.not.ALLOCATED(cell_area)) ALLOCATE(cell_area(ngrid))
     182   IF (.not.ALLOCATED(dx)) ALLOCATE(dx(ngrid))
     183   IF (.not.ALLOCATED(dy)) ALLOCATE(dy(ngrid))
    169184   longitude(:) = plon(:)
    170185   latitude(:) = plat(:)
    … …  
    184199            M_ALBEDO,CST_AL,&
    185200            M_TSURF,M_EMISS,M_CO2ICE,&
    186             M_GW,M_Z0,CST_Z0,&
     201            M_GW,M_Z0,&
    187202            M_H2OICE,&
    188203            phisfi_val)
    … …  
    195210   INTEGER, INTENT(IN) :: ips,ipe,jps,jpe,JULYR,MARS_MODE
    196211   INTEGER :: i,j,subs,nlast
    197    REAL, INTENT(IN  ) :: CST_AL, phisfi_val, CST_Z0
     212   REAL, INTENT(IN  ) :: CST_AL, phisfi_val
    198213   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: &
    199214     M_ALBEDO,M_TSURF,M_EMISS,M_CO2ICE,M_H2OICE,M_Z0
    … …  
    334349            M_Q2,M_WSTAR)
    335350
    336    !use turb_mod, only: q2,wstar,turb_resolved
     351   use turb_mod, only: q2,wstar,turb_resolved
    337352
    338353   implicit none
    … …  
    340355   INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme
    341356   INTEGER, INTENT(IN) :: ips,ipe,jps,jpe
    342    INTEGER :: i,j,subs,nlast     
     357   INTEGER :: i,j,subs,nlast
    343358   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: M_WSTAR
    344359   REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(IN) :: M_Q2
    345360   LOGICAL, INTENT(IN ) :: RESTART,isles
    346361
    347    !! Nothing is done
     362   turb_resolved =.true.
     363   print*,'turb_resolved :',isles
     364
     365   DO j = jps,jpe
     366   DO i = ips,ipe
     367
     368     !-----------------------------------!
     369     ! 1D subscript for physics "cursor" !
     370     !-----------------------------------!
     371     subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
     372
     373     !PBL wind variance
     374     IF (.not. restart) THEN
     375      q2(subs,:) = 1.E-6
     376      wstar(subs)=0.
     377     ELSE
     378      q2(subs,:)=M_Q2(i,:,j)
     379      wstar(subs)=M_WSTAR(i,j)
     380     ENDIF
     381
     382   ENDDO
     383   ENDDO
     384
     385   !!---------------------!!
     386   !! OUTPUT FOR CHECKING !!
     387   !!---------------------!!
     388   nlast = (ipe-ips+1)*(jpe-jps+1)
     389   print*,"check: q2",q2(1,1),q2(nlast,kme+1)
     390   print*,"check: wstar",wstar(1),wstar(nlast)
    348391
    349392END SUBROUTINE update_inputs_physiq_turb

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_venus_lmd_new/update_outputs_physiq_mod.F

    r1596 r1724  
    118118SUBROUTINE update_outputs_physiq_turb( &
    119119            ims,ime,jms,jme,kms,kme,&
    120             ips,ipe,jps,jpe,&
     120            ips,ipe,jps,jpe,kps,kpe,&
    121121            M_Q2,M_WSTAR,&
    122122            HFMAX,ZMAX,USTM,HFX)
    123123
    124    !use turb_mod, only: q2,wstar,ustar,sensibFlux,&
     124   use turb_mod, only: q2,wstar,yustar,sens!,&
    125125   !                     hfmax_th,zmax_th
    126126
    … …  
    128128
    129129   INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme
    130    INTEGER, INTENT(IN) :: ips,ipe,jps,jpe
    131    INTEGER :: i,j,subs   
     130   INTEGER, INTENT(IN) :: ips,ipe,jps,jpe,kps,kpe
     131   INTEGER :: i,j,subs
    132132   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)  :: &
    133133     M_WSTAR,HFMAX,ZMAX,USTM,HFX
    134134   REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(INOUT) :: M_Q2
    135135
     136   DO j = jps,jpe
     137   DO i = ips,ipe
     138
     139     !-----------------------------------!
     140     ! 1D subscript for physics "cursor" !
     141     !-----------------------------------!
     142     subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
     143
     144     !-------------------------------------------------------!
     145     ! Save key variables for restart and output and nesting ! 
     146     !-------------------------------------------------------!
     147     M_Q2(i,kps:kpe+1,j) = q2(subs,:)
     148     M_WSTAR(i,j) = wstar(subs)
     149
     150     !! output only (arrays already in phys modules)
     151     !HFMAX(i,j) = HFMAX_TH(subs)
     152     !ZMAX(i,j) = ZMAX_TH(subs)
     153     USTM(i,j) = yustar(subs)
     154     HFX(i,j) = sens(subs)
     155
     156  ENDDO
     157  ENDDO
    136158
    137159END SUBROUTINE update_outputs_physiq_turb
    … …  
    144166           SWDOWNZ,TAU_DUST,QSURFDUST,&
    145167           MTOT,ICETOT,TAU_ICE,&
    146            HR_SW,HR_LW,&
     168           HR_SW,HR_LW,HR_DYN,DT,DTRAD,DTVDF,DTAJS,&
    147169           RDUST,VMR_ICE,RICE)
    148170
    … …  
    155177     MTOT,ICETOT,TAU_ICE
    156178  REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: &
    157      HR_SW,HR_LW,RDUST,VMR_ICE,RICE
     179     HR_SW,HR_LW,HR_DYN,RDUST,VMR_ICE,RICE,DT,DTRAD,DTVDF,DTAJS
    158180  INTEGER :: i,j,subs
    159181 
    … …  
    171193   HR_SW(i,kps:kpe,j) = comm_HR_SW(subs,kps:kpe)
    172194   HR_LW(i,kps:kpe,j) = comm_HR_LW(subs,kps:kpe)
     195   HR_DYN(i,kps:kpe,j) = comm_HR_DYN(subs,kps:kpe)
     196   DT(i,kps:kpe,j) = comm_DT(subs,kps:kpe)
     197   DTRAD(i,kps:kpe,j) = comm_DT_RAD(subs,kps:kpe)
     198   DTVDF(i,kps:kpe,j) = comm_DT_VDF(subs,kps:kpe)
     199   DTAJS(i,kps:kpe,j) = comm_DT_AJS(subs,kps:kpe)
    173200
    174201   ENDDO

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/dynphy_wrf_void_lmd_new/update_outputs_physiq_mod.F

    r1590 r1724  
    7171            SWDOWNZ,TAU_DUST,QSURFDUST,&
    7272            MTOT,ICETOT,TAU_ICE,&
    73             HR_SW,HR_LW,&
     73            HR_SW,HR_LW,HR_DYN,DT,DTRAD,DTVDF,DTAJS,&
    7474            RDUST,VMR_ICE,RICE)
    7575
    … …  
    8080     MTOT,ICETOT,TAU_ICE
    8181   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: &
    82      HR_SW,HR_LW,RDUST,VMR_ICE,RICE
     82     HR_SW,HR_LW,RDUST,VMR_ICE,RICE,HR_DYN,DT,DTRAD,DTVDF,DTAJS
    8383   INTEGER :: i,j,subs
    8484

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F.new

    r1636 r1724  
    4040        isfflx, diff_opt, km_opt, &
    4141        HISTORY_INTERVAL, &
    42         HR_SW,HR_LW,SWDOWNZ,&
     42        HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,DT_AJS,SWDOWNZ,&
    4343        TAU_DUST,RDUST,QSURFDUST,&
    4444        MTOT,ICETOT,VMR_ICE,TAU_ICE,RICE,&
    … …  
    109109REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
    110110     RTHBLTEN,RUBLTEN,RVBLTEN, &
    111      HR_SW,HR_LW,RDUST,VMR_ICE,RICE
     111     HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,DT_AJS,RDUST,VMR_ICE,RICE
    112112REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(INOUT ) :: &
    113113     M_Q2
    … …  
    138138   ! ------> inputs:
    139139   INTEGER :: ngrid,nlayer,nq,nsoil
    140    REAL :: pday,ptime,MY 
     140   REAL*8 :: pday,ptime
     141   REAL :: MY 
    141142   REAL :: phisfi_val
    142143   LOGICAL :: firstcall,lastcall
    143144   ! ----------
    144    REAL,DIMENSION(:,:),ALLOCATABLE :: pplev,pplay,pphi,pu,pv,pt,flxw
    145    REAL,DIMENSION(:,:,:),ALLOCATABLE :: pq 
     145   REAL*8,DIMENSION(:,:),ALLOCATABLE :: pplev,pplay,pphi,pu,pv,pt,flxw
     146   REAL*8,DIMENSION(:,:,:),ALLOCATABLE :: pq 
    146147
    147148   ! <------ outputs:
    148149   !     physical tendencies
    149    REAL,DIMENSION(:),ALLOCATABLE :: pdpsrf
    150    REAL,DIMENSION(:,:),ALLOCATABLE :: pdu,pdv,pdt
    151    REAL,DIMENSION(:,:,:),ALLOCATABLE :: pdq
     150   REAL*8,DIMENSION(:),ALLOCATABLE :: pdpsrf
     151   REAL*8,DIMENSION(:,:),ALLOCATABLE :: pdu,pdv,pdt,pdtheta
     152   REAL*8,DIMENSION(:,:,:),ALLOCATABLE :: pdq
    152153   ! ... intermediate arrays
    153154   REAL, DIMENSION(:), ALLOCATABLE  :: &
    … …  
    158159   ! Additional control variables
    159160   INTEGER :: sponge_top,relax,ips,ipe,jps,jpe,kps,kpe
    160    REAL :: elaps, ptimestep
     161   REAL :: elaps
     162   REAL*8 :: ptimestep
    161163   INTEGER :: test
    162164   REAL :: wappel_phys
    … …  
    173175             dp_save
    174176      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
    175              du_save, dv_save, dt_save
     177             du_save, dv_save, dt_save,dtheta_save
    176178      REAL, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: &
    177179             dq_save     
    … …  
    180182             dp_save
    181183      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
    182              du_save, dv_save, dt_save
     184             du_save, dv_save, dt_save,dtheta_save
    183185      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
    184186             dq_save     
    … …  
    191193
    192194   !! arguments to physiq
    193    REAL,ALLOCATABLE :: zpk_omp(:,:)
    194    REAL,ALLOCATABLE :: zphis_omp(:) ! surface geopotential
    195    REAL,ALLOCATABLE :: presnivs_omp(:) ! approximate pressure of atm. layers 
    196    REAL,ALLOCATABLE :: zrfi_omp(:,:) ! relative wind vorticity, in s-1
    197 
     195   CHARACTER(len=20),ALLOCATABLE :: tname(:) ! tracer names
     196   REAL*8,ALLOCATABLE :: zpk_omp(:,:)
     197   REAL*8,ALLOCATABLE :: zphis_omp(:) ! surface geopotential
     198   REAL*8,ALLOCATABLE :: presnivs_omp(:) ! approximate pressure of atm. layers 
     199   REAL*8,ALLOCATABLE :: zrfi_omp(:,:) ! relative wind vorticity, in s-1
     200   REAL :: tk1,tk2
    198201!==================================================================
    199202! CODE
    … …  
    301304   ALLOCATE(dv_save(ngrid,nlayer,max_dom))
    302305   ALLOCATE(dt_save(ngrid,nlayer,max_dom))
     306   ALLOCATE(dtheta_save(ngrid,nlayer,max_dom))
    303307   ALLOCATE(dq_save(ngrid,nlayer,nq,max_dom))
    304308   dp_save(:,:)=0.    !! initialize these arrays ...
    … …  
    306310   dv_save(:,:,:)=0.
    307311   dt_save(:,:,:)=0.
     312   dtheta_save(:,:,:)=0.
    308313   dq_save(:,:,:,:)=0.
    309314ENDIF
    … …  
    319324ALLOCATE(dv_save(ngrid,nlayer))
    320325ALLOCATE(dt_save(ngrid,nlayer))
     326ALLOCATE(dtheta_save(ngrid,nlayer))
    321327ALLOCATE(dq_save(ngrid,nlayer,nq))
    322328ENDIF
    … …  
    325331dv_save(:,:)=0.
    326332dt_save(:,:)=0.
     333dtheta_save(:,:)=0.
    327334dq_save(:,:,:)=0.
    328335flag_first_restart=.false.
    … …  
    369376! ALLOCATE !
    370377!----------!
     378IF (.not.ALLOCATED(tname)) ALLOCATE(tname(nq))
    371379!-------------------------------------------------------------------------------!
    372380! outputs:                                                                      !       
    … …  
    381389ALLOCATE(pdv(ngrid,nlayer))
    382390ALLOCATE(pdt(ngrid,nlayer))
     391ALLOCATE(pdtheta(ngrid,nlayer))
    383392ALLOCATE(pdq(ngrid,nlayer,nq))
    384393!!!
    … …  
    392401pdv(:,:)=dv_save(:,:,id)
    393402pdt(:,:)=dt_save(:,:,id)
     403pdtheta(:,:)=dtheta_save(:,:,id)
    394404pdq(:,:,:)=dq_save(:,:,:,id)
    395405#else
     406print*,'else'
    396407pdpsrf(:)=dp_save(:)
    397408pdu(:,:)=du_save(:,:)
    398409pdv(:,:)=dv_save(:,:)
    399410pdt(:,:)=dt_save(:,:)
     411pdtheta(:,:)=dtheta_save(:,:)
    400412pdq(:,:,:)=dq_save(:,:,:)
    401413#endif
    … …  
    442454  !! tracers' name
    443455  PRINT *,'** ',planet_type,'** TRACERS NAMES'
    444   CALL update_inputs_physiq_tracers(nq,MARS_MODE)
     456  CALL update_inputs_physiq_tracers(nq,MARS_MODE,tname)
    445457  !! PHYSICS VARIABLES (cf. iniphysiq in LMD GCM)
    446458  !! parameters are defined in the module_model_constants.F WRF routine
    … …  
    504516  ENDIF
    505517ENDIF
    506 
    507 IF (MARS_MODE .EQ. 32) THEN
    508   IF (firstcall .EQV. .true. .and. (.not. restart)) THEN
    509       q_prof(:,7) = 0.95
    510 !! traceurs(7)   = 'co2'
    511   ENDIF
    512 ENDIF
    513 
    514518
    515519IF (firstcall .EQV. .true.) THEN
    … …  
    621625            M_ALBEDO,CST_AL,&
    622626            M_TSURF,M_EMISS,M_CO2ICE,&
    623             M_GW,M_Z0,CST_Z0,&
     627            M_GW,M_Z0,&
    624628            M_H2OICE,&
    625629            phisfi_val)
    … …  
    653657pdv(:,:)=0.
    654658pdt(:,:)=0.
     659pdtheta(:,:)=0.
    655660pdq(:,:,:)=0.
    656661print *, '** ',planet_type,'** CALL TO LMD PHYSICS'
    … …  
    662667            ptime,pday,MY)
    663668!!!
    664 CALL call_physiq(planet_type,ngrid,nlayer,nq,      &
     669CALL call_physiq(planet_type,ngrid,nlayer,nq,tname, &
    665670                       firstcall,lastcall,          &
    666671                       pday,ptime,ptimestep,        &
    … …  
    679684ENDIF
    680685#endif
     686
     687IF (planet_type .eq. "venus" ) THEN
     688  DO j=jps,jpe
     689  DO i=ips,ipe
     690    do k=kps,kpe
     691      subs=(j-jps)*(ipe-ips+1)+(i-ips+1)
     692      tk1=(pt(subs,k)**nu + nu*TT00**nu*log((p1000mb/pplay(subs,k))**rcp))**(1/nu)
     693      tk2=((pt(subs,k) + pdt(subs,k))**nu + nu*TT00**nu*log((p1000mb/pplay(subs,k))**rcp))**(1/nu)
     694      pdtheta(subs,k)=tk2-tk1
     695    enddo
     696  ENDDO
     697  ENDDO
     698ENDIF
    681699
    682700print *, '** ',planet_type,'** CALL TO LMD PHYSICS DONE'
    … …  
    694712DEALLOCATE(zrfi_omp)
    695713
    696 
    697714!---------------------------------------------------------------------------------!
    698715! PHYSIQ TENDENCIES ARE SAVED TO BE SPLIT WITHIN INTERMEDIATE DYNAMICAL TIMESTEPS !
    … …  
    703720dv_save(:,:,id)=pdv(:,:)
    704721dt_save(:,:,id)=pdt(:,:)
     722dtheta_save(:,:,id)=pdtheta(:,:)
    705723dq_save(:,:,:,id)=pdq(:,:,:)
    706724#else
    … …  
    709727dv_save(:,:)=pdv(:,:)
    710728dt_save(:,:)=pdt(:,:)
     729dtheta_save(:,:)=pdtheta(:,:)
    711730dq_save(:,:,:)=pdq(:,:,:)
    712731#endif
    … …  
    737756CALL update_outputs_physiq_turb( &
    738757            ims,ime,jms,jme,kms,kme,&
    739             ips,ipe,jps,jpe,&
     758            ips,ipe,jps,jpe,kps,kpe,&
    740759            M_Q2,M_WSTAR,&
    741760            HFMAX,ZMAX,USTM,HFX)
    … …  
    746765            SWDOWNZ,TAU_DUST,QSURFDUST,&
    747766            MTOT,ICETOT,TAU_ICE,&
    748             HR_SW,HR_LW,&
     767            HR_SW,HR_LW,HR_DYN,DDT,DT_RAD,DT_VDF,DT_AJS,&
    749768            RDUST,VMR_ICE,RICE)
    750769!!!
     770print*,"update_outputs_physiq_diag"
     771
     772
    751773ENDIF call_physics
    752774
    … …  
    770792! --is the one calculated during the last call to physics          !
    771793!------------------------------------------------------------------!
    772 
     794!print*,'pdt',pdt(1,1),pdt(1,nlayer)
     795!print*,'exner',exner(1,:,1)
    773796DO j = jps,jpe
    774797DO i = ips,ipe
    … …  
    778801    ! zonal wind
    779802  RUBLTEN(i,kps:kpe,j) = pdu(subs,kps:kpe)
    780 
    781803    ! meridional wind
    782804  RVBLTEN(i,kps:kpe,j) = pdv(subs,kps:kpe)
    783 
    784805    ! potential temperature
    785806    ! (dT = dtheta * exner for isobaric coordinates or if pressure variations are negligible)
    786   RTHBLTEN(i,kps:kpe,j) = pdt(subs,kps:kpe) / exner(i,kps:kpe,j)
    787 
     807  IF (planet_type .eq. "venus" ) THEN
     808    RTHBLTEN(i,kps:kpe,j) = pdtheta(subs,kps:kpe)
     809  ELSE
     810   RTHBLTEN(i,kps:kpe,j) = pdt(subs,kps:kpe) / exner(i,kps:kpe,j)
     811  ENDIF
    788812    ! update surface pressure (cf CO2 cycle in physics)
    789813    ! here dt is needed
    790814  PSFC(i,j)=PSFC(i,j)+pdpsrf(subs)*dt
    791 
    792815    ! tracers
    793816  SCALAR(i,kps:kpe,j,1)=0.
    … …  
    807830ENDDO
    808831ENDDO
    809 
    810832DEALLOCATE(pdpsrf)
    811833DEALLOCATE(pdu)
    … …  
    813835DEALLOCATE(pdt)
    814836DEALLOCATE(pdq)
    815 
     837DEALLOCATE(pdtheta)
    816838!!*****!!
    817839!! END !!

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/share/Makefile

    r11 r1724  
    1717        module_optional_si_input.o      \
    1818        module_compute_geop.o           \
    19         module_soil_pre.o
     19        module_soil_pre.o \
     20        cpdet_mod.o
    2021
    2122OBJS    =               \

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/share/module_model_constants.F

    r1604 r1724  
    88
    99   !  A really small number.
    10 
     10   character(len=15) :: planet
    1111   REAL    , PARAMETER :: epsilon         = 1.E-15
    1212
    … …  
    183183    REAL  :: z_scale        ! scale height
    184184
     185    REAL :: TT00
     186    REAL :: nu
    185187 CONTAINS
    186188  SUBROUTINE init_planet_constants
    187189    IMPLICIT NONE
    188     character(len=15) :: planet
    189190    Call initial_config
    190     planet=model_config_rec%planet
    191191    !write(*,*) "unknown planet type", planet
    192192    !stop
    … …  
    239239     wvolcapa      = 1.e6       ! volumetric capacity of soil (new soil model)
    240240#endif
     241     TT00=0.
     242     nu=0.
    241243!#endif
    242244    ELSE if ( planet == "venus") then
    … …  
    255257     womeg        = 0.2992549E-06
    256258     wdaysec      = 0.1008707E08
     259     TT00=460.
     260     nu=0.35
    257261    ELSE if ( planet == "prescribed") then
    258262     open(17,file='planet_constant',form='formatted',status='old')
    … …  
    268272     p1000mb=p0
    269273     reradius=1./rad
     274     TT00=0.
     275     nu=0.
    270276    Else
    271277     write(*,*) "unknown planet type:", planet
    … …  
    295301    write(*,*) "t0 = ", t0
    296302    write(*,*) "p0 = ", p0
     303    write(*,*) "T0 = ", TT00
     304    write(*,*) "nu = ", nu
    297305  END SUBROUTINE init_planet_constants
    298306 END MODULE module_model_constants
Note: See TracChangeset for help on using the changeset viewer.