Ignore:
Timestamp:
Jun 4, 2015, 4:23:32 PM (9 years ago)
Author:
slebonnois
Message:

SL: update of the Venus GCM, + corrections on routines used for newstart/start2archive for Titan and Venus, + some modifications on tools

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.VENUS/libf/phyvenus/physiq.F

    r1310 r1442  
    6464      USE write_field_phy
    6565      USE iophy
    66       use cpdet_mod, only: cpdet, t2tpot
     66      USE cpdet_mod, only: cpdet, t2tpot
    6767      USE chemparam_mod
    6868      USE conc
    … …  
    240240      EXTERNAL conduction
    241241      EXTERNAL molvis
     242      EXTERNAL moldiff_red
     243
    242244c
    243245c Variables locales
    … …  
    260262      REAL zdtime, zlongi
    261263c
    262       INTEGER i, k, iq, ig, j, ll
     264      INTEGER i, k, iq, ig, j, ll, ilon, ilat, ilev
    263265c
    264266      REAL zphi(klon,klev)
    … …  
    301303      real d_v_molvis(klon,klev)     ! (m/s) /s
    302304
     305c Tendencies due to molecular diffusion
     306      real d_q_moldif(klon,klev,nqmax)
     307
    303308c
    304309c Variables liees a l'ecriture de la bande histoire physique
    … …  
    322327      REAL :: tr_seri(klon,klev,nqmax)
    323328      REAL :: d_tr(klon,klev,nqmax)
     329
     330c Champ de modification de la temperature par rapport a VIRAII
     331      REAL delta_temp(klon,klev)
     332c      SAVE delta_temp
     333      REAL mat_dtemp(33,50)
     334      SAVE mat_dtemp
    324335
    325336c Variables tendance sedimentation
    … …  
    493504C TRACEURS
    494505C source dans couche limite
    495          source = 0.0 ! pas de source, pour l'instant
     506         source(:,:) = 0.0 ! pas de source, pour l'instant
    496507c---------
    497508
    … …  
    510521            rnew(ig,j)=R
    511522            cpnew(ig,j)=cpdet(tmoy(j))
    512 c            print*, ' physique  l503'
    513 c            print*,  j, cpdet(tmoy(j))
    514              mmean(ig,j)=RMD
     523            mmean(ig,j)=RMD
    515524            akknew(ig,j)=1.e-4
    516525          enddo
    … …  
    522531         call compo_hedin83_init2
    523532      ENDIF
    524       if (callnlte) call nlte_setup
    525       if(callnirco2.and.(nircorr.eq.1)) call nir_leedat         
     533      if (callnlte.and.nltemodel.eq.2) call nlte_setup
     534      if (callnirco2.and.nircorr.eq.1) call nir_leedat         
    526535c---------
    527536     
    … …  
    627636      endif
    628637           
    629       if ((nlon .GT. 1) .AND. ok_chem) then
     638      if ((nlon .GT. 1) .AND. (ok_chem.OR.ok_cloud)) then
    630639c !!! DONC 3D !!!
    631640        CALL chemparam_ini()
    … …  
    636645        CALL cloud_ini(nlon,nlev)
    637646      endif
     647
     648c======================================================================
     649c      Lecture du fichier DeltaT
     650c======================================================================
     651
     652c     ATTENTION tout ce qui suit est pour un 48*32*50
     653
     654      if (ok_deltatemp) then
     655
     656      print*,'lecture de VenusDeltaT.txt '
     657      open(99, form = 'formatted', status = 'old', file =
     658     &     'VenusDeltaT.dat')
     659      print*,'Ouverture de VenusDeltaT.txt '
     660
     661      DO ilev = 1, klev
     662         read(99,'(33(1x,e13.6))') (mat_dtemp(ilat,ilev),ilat=1,33)
     663         print*,'lecture de VenusDeltaT.txt ligne:',ilev
     664      ENDDO
     665     
     666      close(99)
     667      print*,'FIN lecture de VenusDeltaT.txt ok.'
     668
     669      DO k = 1, klev
     670      DO i = 1, klon     
     671         ilat=(rlatd(i)/5.625) + 17.
     672         delta_temp(i,k)=mat_dtemp(INT(ilat),k)
     673      ENDDO
     674      ENDDO
     675     
     676      endif
    638677       
    639678      ENDIF ! debut
    640 c====================================================================
     679c======================================================================
    641680c======================================================================
    642681
    … …  
    830869! Case 3: Full chemistry and/or clouds
    831870
    832          call phytrac_chimie(                 
     871         if (ok_deltatemp) then
     872!           PRINT*,'Def de delta_temp'
     873           DO k = 1, klev
     874           DO i = 1, klon     
     875             ilat=(rlatd(i)/5.625) + 17.
     876!             PRINT*,INT(ilat),rlatd(i),mat_dtemp(INT(ilat),k)
     877             delta_temp(i,k)=mat_dtemp(INT(ilat),k)
     878           ENDDO
     879           ENDDO
     880     
     881         endif
     882
     883         if (ok_deltatemp) then
     884! Utilisation du champ de temperature modifie         
     885           call phytrac_chimie(                             
    833886     I             debut,
    834887     I             gmtime,
    835888     I             nqmax,
    836      I             nlon,
     889     I             klon,
    837890     I             rlatd,
    838891     I             rlond,
    839892     I             nlev,
    840893     I             dtime,
    841      I             t_seri,pplay,
    842      O             tr_seri,
    843      O             NBRTOT,
    844      O             WH2SO4,
    845      O             rho_droplet)
     894     I             t_seri+delta_temp,
     895     I             pplay,
     896     O             tr_seri)
     897         else
     898         
     899           call phytrac_chimie(                             
     900     I             debut,
     901     I             gmtime,
     902     I             nqmax,
     903     I             klon,
     904     I             rlatd,
     905     I             rlond,
     906     I             nlev,
     907     I             dtime,
     908     I             t_seri,
     909     I             pplay,
     910     O             tr_seri)
     911         endif
    846912
    847913c        CALL WriteField_phy('Pression',pplay,nlev)
    … …  
    856922         if (ok_sedim) then
    857923
     924         if (ok_deltatemp) then
     925! Utilisation du champ de temperature modifie 
    858926           CALL new_cloud_sedim(
    859      I               klon,
    860      I               nlev,
    861      I               dtime,
    862      I               pplay,
    863      I               paprs,
    864      I               t_seri,
    865      I               WH2SO4,
    866      I               tr_seri,
    867      I               nqmax,
    868      I               NBRTOT,
    869      I               rho_droplet,
    870      O               Fsedim,
    871      O               d_tr_sed,
    872      O               d_tr_ssed)
    873 
    874           DO k = 1, klev
    875            DO i = 1, klon
    876      
    877 c        WRITE(88,"(11(e15.8,','))") pplay(5,25),
    878 c     &  t_seri(5,25),tr_seri(5,25,i_h2oliq),
    879 c     &  tr_seri(5,25,i_h2o),tr_seri(5,25,i_h2so4liq),
    880 c     &  tr_seri(5,25,i_h2so4),NBRTOT(5,25),WH2SO4(5,25),
    881 c     &  Fsedim(5,25),d_tr_sed(5,25,1),d_tr_sed(5,25,2)
     927     I                 klon,
     928     I               nlev,
     929     I               dtime,
     930     I                 pplay,
     931     I               paprs,
     932     I               t_seri+delta_temp,
     933     I                 tr_seri,
     934     O               d_tr_sed,
     935     O               d_tr_ssed,
     936     I               nqmax,
     937     O                 Fsedim)
     938         else
     939         
     940           CALL new_cloud_sedim(
     941     I                 klon,
     942     I               nlev,
     943     I               dtime,
     944     I                 pplay,
     945     I               paprs,
     946     I               t_seri,
     947     I                 tr_seri,
     948     O               d_tr_sed,
     949     O               d_tr_ssed,
     950     I               nqmax,
     951     O                 Fsedim)
     952
     953         endif
     954
     955        DO k = 1, klev
     956         DO i = 1, klon
    882957
    883958c--------------------
    … …  
    888963        PRINT*,'d_tr_sed Nan?',d_tr_sed(i,k,:),'Temp',t_seri(i,k)
    889964        PRINT*,'lat-lon',i,'level',k,'dtime',dtime
    890         PRINT*,'NBRTOT',NBRTOT(i,k),'F_sed',Fsedim(i,k)
     965        PRINT*,'F_sed',Fsedim(i,k)
    891966        PRINT*,'==============================================='
    892967                d_tr_sed(i,k,:)=0.
    … …  
    901976        Fsedim(i,k)     = Fsedim(i,k) / dtime
    902977     
    903            ENDDO
    904978          ENDDO
     979         ENDDO
    905980     
    906981        Fsedim(:,klev+1) = 0.
    … …  
    9881063         ENDDO
    9891064         ENDDO
     1065   
    9901066         DO iq=1, nqmax
     1067c  AS: changement
     1068c  Pourquoi d_tr_vdf(1,1,iq) et tr_seri(1,1,iq)
     1069c  et pas d_tr_vdf(:,:,iq) tr_seri(:,:,iq)
     1070c  Je vois pas en quoi cltrac ne prendrait en compte que le traceur à la surface et au point 1 en klon
     1071c
     1072
     1073c  Je garde le source(:,iq) parce que je comprend pas sinon source
     1074c   dimension(klon,nqmax) et flux dans cltrac (klon) ???
     1075
     1076c             CALL cltrac(dtime,ycoefh,t_seri,
     1077c     s               tr_seri(1,1,iq),source(:,iq),
     1078c     e               paprs, pplay,delp,
     1079c     s               d_tr_vdf(1,1,iq))
     1080     
    9911081             CALL cltrac(dtime,ycoefh,t_seri,
    992      s               tr_seri(1,1,iq),source,
     1082     s               tr_seri(:,:,iq),source(:,iq),
    9931083     e               paprs, pplay,delp,
    994      s               d_tr_vdf(1,1,iq))
     1084     s               d_tr_vdf(:,:,iq))
     1085     
    9951086             tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
    9961087             d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime          ! /s
     1088
     1089         DO k = 1, klev
     1090         DO i = 1, klon
     1091            tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr_vdf(i,k,iq)
     1092            d_tr_vdf(i,k,iq)= d_tr_vdf(i,k,iq)/dtime          ! /s
    9971093         ENDDO
    998       endif
     1094         ENDDO
     1095         
     1096         ENDDO !nqmax
     1097
     1098       endif
    9991099
    10001100      IF (if_ebil.ge.2) THEN
    … …  
    10841184           d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime  ! /s
    10851185         endif
    1086 
    10871186      endif
    10881187
    … …  
    11041203      END IF
    11051204
     1205
    11061206c====================================================================
    11071207c RAYONNEMENT
    11081208c====================================================================
    11091209
    1110 c-----------------------------------------------------------------------
     1210c------------------------------------
    11111211c    . Compute radiative tendencies :
    11121212c------------------------------------
    … …  
    11181218c     PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
    11191219           
    1120 c Calcul pour Cp rnew et mmean avec traceurs (Cp independant de T !! )
    1121        IF(callnlte.or.callthermos) THEN                                 
     1220
     1221c------------------------------------
     1222c    . Compute mean mass, cp and R :
     1223c------------------------------------
     1224
     1225      if(callthermos) then
     1226         call concentrations2(pplay,t_seri,d_t,tr_seri, nqmax,
     1227     &                        pdtphys)
     1228
     1229      endif
     1230
     1231
     1232cc!!! ADD key callhedin
     1233
     1234      IF(callnlte.or.callthermos) THEN                                 
    11221235         call compo_hedin83_mod(pplay,rmu0,   
    11231236     &           co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
     1237
     1238         IF(ok_chem) then
     1239 
     1240CC  !! GG : Using only mayor species tracers abundances to compute NLTE heating/cooling
     1241
     1242CC               Conversion [mmr] ---> [vmr]
     1243       
     1244                 co2vmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_co2)*
     1245     &                             mmean(1:nlon,1:nlev)/M_tr(i_co2)
     1246                 covmr_gcm(:,:)  = tr_seri(1:nlon,1:nlev,i_co)*
     1247     &                              mmean(1:nlon,1:nlev)/M_tr(i_co)
     1248                 ovmr_gcm(:,:)   = tr_seri(1:nlon,1:nlev,i_o)*
     1249     &                             mmean(1:nlon,1:nlev)/M_tr(i_o)
     1250                 n2vmr_gcm(:,:)   = tr_seri(1:nlon,1:nlev,i_n2)*
     1251     &                             mmean(1:nlon,1:nlev)/M_tr(i_n2)
     1252
     1253         ENDIF
     1254
    11241255       ENDIF   
    11251256
    1126       if(callthermos) then
    1127          call concentrations2(pplay,t_seri,d_t,co2vmr_gcm, n2vmr_gcm,
    1128      &          covmr_gcm, ovmr_gcm,nvmr_gcm,pdtphys)
    1129       endif
    1130 
    1131 
    1132 c          NLTE cooling from CO2 emission
    1133 c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     1257c
     1258c   NLTE cooling from CO2 emission
     1259c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    11341260
    11351261        IF(callnlte) THEN                                 
    11361262            if(nltemodel.eq.0.or.nltemodel.eq.1) then
    1137                 CALL nltecool(klon, klev, pplay*9.869e-6, t_seri,
    1138      $                    co2vmr_gcm,n2vmr_gcm, covmr_gcm, ovmr_gcm,
    1139      $                    d_t_nlte)
     1263                CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
     1264     $                    tr_seri, d_t_nlte)
    11401265            else if(nltemodel.eq.2) then                               
    1141                 CALL nlte_tcool(klon,klev,pplay*9.869e-6,             
     1266               CALL nlte_tcool(klon,klev,pplay*9.869e-6,             
    11421267     $               t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
    11431268     $               ovmr_gcm,d_t_nlte,ierr_nlte,varerr )
    … …  
    11621287     $        CALL nlthermeq(klon, klev, paprs, pplay)
    11631288
    1164 
    1165 c          LTE radiative transfert / solar / IR matrix
    1166 c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1167 
     1289c
     1290c       LTE radiative transfert / solar / IR matrix
     1291c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    11681292      CALL radlwsw
    11691293     e            (dist, rmu0, fract, zzlev,
    … …  
    11711295
    11721296
    1173 c          CO2 near infrared absorption
    1174 c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     1297c       CO2 near infrared absorption
     1298c      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    11751299
    11761300        d_t_nirco2(:,:)=0.
    11771301        if (callnirco2) then
    1178            call nirco2abs (klon, klev, pplay, dist, nqmax, qx,
    1179      .                 rmu0, fract, d_t_nirco2,
    1180      .                 co2vmr_gcm, ovmr_gcm)
     1302           call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
     1303     .                 rmu0, fract, d_t_nirco2)
    11811304        endif
    11821305
    … …  
    12091332        IF (callthermos) THEN
    12101333
    1211 c        call thermosphere(zplev,zplay,dist_sol,
    1212 c     $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
    1213 c     &     pt,pq,pu,pv,pdt,pdq,
    1214 c     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
    1215 
    1216            call euvheat(klon, klev, t_seri,paprs,pplay,zzlay,
    1217      $          rmu0,pdtphys,gmtime,rjourvrai,
    1218 C     $                 pq,pdq,zdteuv)
    1219      $          co2vmr_gcm, n2vmr_gcm, covmr_gcm,
    1220      $          ovmr_gcm,nvmr_gcm,d_t_euv )
    1221 
     1334c           call euvheat(klon, klev,t_seri,paprs,pplay,zzlay,
     1335c     $          rmu0,pdtphys,gmtime,rjourvrai, co2vmr_gcm, n2vmr_gcm,
     1336c     $          covmr_gcm, ovmr_gcm,d_t_euv )
     1337           call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
     1338     $         rmu0,pdtphys,gmtime,rjourvrai,
     1339     $         tr_seri, d_tr, d_t_euv )
     1340               
    12221341           DO k=1,klev
    12231342              DO ig=1,klon
    … …  
    12621381            call molvis(klon, klev, pdtphys,
    12631382     $            pplay,paprs,t_seri,
    1264      $            v,tsurf,zzlev,zzlay,d_u_molvis)
     1383     $            v,tsurf,zzlev,zzlay,d_v_molvis)
    12651384
    12661385            DO k=1,klev
    … …  
    12711390               ENDDO
    12721391            ENDDO
    1273 
    1274          ENDIF  ! callthermos
     1392        ENDIF
     1393
     1394
     1395!  --  MOLECULAR DIFFUSION ---
     1396
     1397          d_q_moldif(:,:,:)=0
     1398
     1399         IF (callthermos .and. ok_chem) THEN
     1400
     1401             call moldiff_red(klon, klev, nqmax,
     1402     &                   pplay,paprs,t_seri, tr_seri, pdtphys,
     1403     &                   zzlay,d_t_euv,d_t_conduc,d_q_moldif)
     1404
     1405
     1406! --- update tendencies tracers ---
     1407
     1408          DO iq = 1, nqmax
     1409           DO k=1,klev
     1410              DO ig=1,klon
     1411                tr_seri(ig,k,iq)= tr_seri(ig,k,iq)+
     1412     &                           d_q_moldif(ig,k,iq)*dtime ! [Kg/kg]?
     1413              ENDDO
     1414            ENDDO
     1415           ENDDO
     1416           
     1417
     1418         ENDIF  ! callthermos & ok_chem
    12751419
    12761420c====================================================================
Note: See TracChangeset for help on using the changeset viewer.