Changeset 1301 for trunk/LMDZ.VENUS/libf


Ignore:
Timestamp:
Jun 26, 2014, 12:25:31 PM (11 years ago)
Author:
slebonnois
Message:

SL: many bug corrections in phyvenus, some cleaning, and a new ksi matrix format for Venus IR

Location:
trunk/LMDZ.VENUS/libf/phyvenus
Files:
23 edited

Legend:

Unmodified
Added
Removed

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

    r1017 r1301  
    147147              ENDIF
    148148c ----------------------------
    149 c TEST --- PAS DE MELANGE DE U
     149c TEST --- PAS DE MELANGE DE U ni Q
    150150c             zalpha=0.
    151151c ----------------------------

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

    r1017 r1301  
    368368         ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
    369369         ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
    370 c --
    371 c VENUS: on met le coefficient K =0 au-dessus de 40 km (2.5e5 Pa)
    372          do i=1,knon
    373          do k=2,klev
    374            if (ypaprs(i,k).lt.2.5e5) then
    375              ycoefm(i,k)=1.e-7
    376              ycoefh(i,k)=1.e-7
    377            endif
    378          enddo
    379          enddo
    380 c --
    381370
    382371c-------------------------------------------------

  • trunk/LMDZ.VENUS/libf/phyvenus/comcstVE.h

    r892 r1301  
    22! INCLUDE comcstVE.h
    33
    4         integer nnuve,nbztopve,nbpsve,nbmat
    5         parameter (nnuve=68)   ! fichiers Vincent et Bullock
    6         parameter (nbpsve=16)   ! number of psurf in Vincent's matrixes
    7         parameter (nbztopve=4) ! number of ztop  in Vincent's matrixes
    8         parameter (nbmat=nbztopve*nbpsve) ! number of matrixes in Vincent's file
     4        integer nnuve,nbmat
     5        parameter (nnuve=68)  ! fichiers Vincent et Bullock
     6        parameter (nbmat=210) ! Max number of matrixes in Vincent's file
    97
    10       common/comcstVE/al,bl,ztopve,psurfve
     8      common/comcstVE/al,bl,nlatve,indexve,nbpsve,nbszave,               &
     9     & psurfve,szave
    1110
    12       real   al(nnuve),bl(nnuve)      ! for Planck luminances calculations
    13       real   ztopve(nbmat)   ! cloud top altitude in matrixes (km)
    14       real   psurfve(nbmat)            ! surface pressure in matrixes (Pa)
     11      real   al(nnuve),bl(nnuve)     ! for Planck luminances calculations
     12! Structure of the ksi matrixes
     13      integer nlatve,indexve(5),nbpsve(5),nbszave(5)
     14      real   psurfve(16,5)           ! surface pressure in matrixes (Pa)
     15      real   szave(8,5)              ! converted in mu0
    1516
    1617

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

    r101 r1301  
    3030c paprs---input-R-Pressure in semi layers    (Pa)
    3131c pplay---input-R-Pressure model-layers      (Pa)
    32 c pgeop---input-R-Geopotential model layers      (m)
     32c pgeop---input-R-Geopotential model layers (reference to ground)
    3333c pn2-----input-R-Brunt-Vaisala freq.^2 at 1/2 layers
    3434c t-------input-R-temperature (K)
    … …  
    6666c ================
    6767c
    68 c zgeom-----R: Altitude of layer above ground
     68c zgeom-----R: Altitude (m) of layer above ground (from top to bottom)
    6969c pt, pu, pv --R: t u v from top to bottom
    7070c pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom)
    … …  
    137137      DO k = klev, 1, -1
    138138      DO i = 1, klon
    139          zgeom(i,k) = pgeop(i,klev-k+1)
     139         zgeom(i,k) = pgeop(i,klev-k+1)/RG
    140140         zn2(i,k)   = pn2(i,klev-k+1)
    141141      ENDDO

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

    r1048 r1301  
    9393      REAL zpic(imar+1,jmar),zval(imar+1,jmar)
    9494      real num_tot(2200,1100),num_lan(2200,1100)
    95 c
     95
    9696      REAL a(2200),b(2200),c(1100),d(1100)
    97 c
     97
     98c pas defini puisque pas de physique dans newstart...
     99      RPI=2.*ASIN(1.)
     100      RA=6051300.
     101
    98102      print *,' parametres de l orographie a l echelle sous maille'
    99103

  • trunk/LMDZ.VENUS/libf/phyvenus/ini_histmth.h

    r902 r1301  
    151151     .                32, "ave(X)", zsto,zout)
    152152c
    153          CALL histdef(nid_mth, "fluxvdf", "PBL net flux","W/m2",
    154      .                iim,jj_nb,nhori, klev,1,klev,nvert,
    155      .                32, "ave(X)", zsto,zout)
    156 c
    157          CALL histdef(nid_mth, "fluxdyn", "Dyn. net flux","W/m2",
    158      .                iim,jj_nb,nhori, klev,1,klev,nvert,
    159      .                32, "ave(X)", zsto,zout)
    160 c
    161          CALL histdef(nid_mth, "fluxajs", "Dry adj. net flux","W/m2",
    162      .                iim,jj_nb,nhori, klev,1,klev,nvert,
    163      .                32, "ave(X)", zsto,zout)
     153c        CALL histdef(nid_mth, "fluxvdf", "PBL net flux","W/m2",
     154c    .                iim,jj_nb,nhori, klev,1,klev,nvert,
     155c    .                32, "ave(X)", zsto,zout)
     156c
     157c        CALL histdef(nid_mth, "fluxdyn", "Dyn. net flux","W/m2",
     158c    .                iim,jj_nb,nhori, klev,1,klev,nvert,
     159c    .                32, "ave(X)", zsto,zout)
     160c
     161c        CALL histdef(nid_mth, "fluxajs", "Dry adj. net flux","W/m2",
     162c    .                iim,jj_nb,nhori, klev,1,klev,nvert,
     163c    .                32, "ave(X)", zsto,zout)
    164164c
    165165c        CALL histdef(nid_mth, "fluxec", "Cin. net flux","W/m2",

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

    r892 r1301  
    2323C     MODIFICATIONS.
    2424C     --------------
    25 C        version multimatrice (topographie, sommet nuages): 20/12/2006
     25C       
     26c   New ksi matrix: possibility of different cloud model fct of lat   05/2014
    2627C     ------------------------------------------------------------------
    2728C
    … …  
    3031c inputs
    3132      real   psurf(klon)           ! Surface pressure
    32       real   ztop(klon)            ! Altitude of the top of cloud deck (km)
    3333c outputs
    3434      real   ksive(0:klev+1,0:klev+1,nnuve,nbmat)  ! ksi matrixes in Vincent's file
    3535
    3636c local variables
    37       integer i,j,ig,band,pve,nlve
    38       integer mat,Nb,m,Nmat,nl_init,mat0
    39       parameter(nl_init=8)
     37      integer i,j,isza,ips,band,pve,sve,nlve
     38      integer lat,Nb,m,mat
    4039      character*9 tmp1
    4140      character*100 file
     41      CHARACTER*2 str2
    4242      real   lambda(nnuve)            ! wavelenght in table (mu->m, middle of interval)
    4343      real   lambdamin(nnuve),lambdamax(nnuve) ! in microns
    44       real   dlambda                  ! cm-1
     44      real   dlambda                  ! in microns
    4545
    4646      nlve = klev
    … …  
    5353      open(10,file=file)
    5454     
    55       do i=1,nl_init-1
    56          read(10,*)
    57       enddo
    58       read(10,*) (tmp1,i=1,4),Nmat
     55      read(10,*)
     56      read(10,*) nlatve
     57      read(10,*)
    5958
    60       if (nbmat.ne.Nmat) then
    61          write(*,*) 'This is subroutine load_ksi'
    62          print*,'Probleme de dimension entre ksi.txt et le param nbmat'
    63          print*,'Nb matrices = ',nbmat,Nmat
    64          stop
    65       endif
    66 
    67       do mat=1,nbmat
     59      write(*,*) 'This is subroutine load_ksi'
     60      write(*,*) 'Nb of lat bands:',nlatve
     61     
     62      do lat=1,nlatve
     63        read(10,*) !line for lat range
     64        read(10,*) indexve(lat)
     65        read(10,*) nbpsve(lat)
    6866        read(10,*)
    69         read(10,*)
     67        read(10,*) nbszave(lat)
     68        read(10,*)
     69       
     70        do isza=1,nbszave(lat)
     71          do ips=1,nbpsve(lat)
     72         
     73        read(10,*) (tmp1,j=1,3),mat    !line for matrix number
    7074        read(10,*) (tmp1,j=1,2),pve
    71         psurfve(mat) = pve*1.e5  ! pve en bar, psurfve en Pa
    72         read(10,*) (tmp1,j=1,7),ztopve(mat)
    73         ztopve(mat) = ztopve(mat)*1.e-3 ! passage en km
     75        psurfve(ips,lat) = pve*1.e5    ! pve in bar, psurfve in Pa
     76        read(10,*) (tmp1,j=1,3),sve
     77        szave(isza,lat) = cos(sve*RPI/180.) ! conversion in mu0
    7478        read(10,*)
    7579        read(10,*) m,Nb
    7680        if (m.ne.nlve) then
    7781         write(*,*) 'This is subroutine load_ksi'
    78          print*,'Probleme de dimension entre ksi.txt et le param nlve'
     82         print*,'Dimension problem between ksi.txt and nlve'
    7983         print*,'N levels = ',m,nlve
    8084         stop
    … …  
    8286        if (Nb.ne.nnuve) then
    8387         write(*,*) 'This is subroutine load_ksi'
    84          print*,'Probleme de dimension entre ksi.txt et le param nnuve'
     88         print*,'Dimension problem between ksi.txt and nnuve'
    8589         print*,'N freq = ',Nb,nnuve
    8690         stop
    8791        endif
    8892c     Now reading ksi matrix index "mat"
     93        write(str2,'(i2.2)') m+2
    8994        do band=1,Nb
    9095         read(10,*) lambdamin(band),lambdamax(band)
    9196         do i=0,m+1
    92             read(10,'(100e17.9)') (ksive(i,j,band,mat),j=0,m+1) ! sr/µm/cm¯¹
     97            read(10,'('//str2//'e17.9)') (ksive(i,j,band,mat),j=0,m+1) ! no unit
    9398         enddo                  ! i
    9499        enddo                     ! band
    95100c       print*,"Matrice ",mat," lue"
    96 c       print*,"   psurf=",psurfve(mat)," bars, Ztop=",ztopve(mat)," km"
    97       enddo  ! mat
     101c       print*,"   psurf=",psurfve(ips,lat)," bars"
     102        if (mat+1.gt.nbmat) then
     103         write(*,*) 'This is subroutine load_ksi'
     104         print*,'Dimension problem between ksi.txt and nbmat'
     105         print*,'(max number of matrixes)'
     106         print*,'nbmat (in comcstVE.h) should be raised'
     107         stop
     108        endif
     109
     110          enddo    ! ips
     111        enddo      ! isza
     112      enddo        ! lat
     113     
     114      write(*,*) 'Total nb of matrixes:',mat
    98115     
    99116      close(10)
    100117
    101 c longueur d'onde centrale et largeur de chaque bande
     118c central wavelength and wavelength width
    102119      do band=1,nnuve
    103          lambda(band)=(lambdamin(band)+lambdamax(band))/2.*1.e-6   ! en m
    104          dlambda     =(1./lambdamin(band)-1./lambdamax(band))*1.e4 ! en cm-1
     120         lambda(band)=(lambdamin(band)+lambdamax(band))/2.*1.e-6   ! in m
     121         dlambda     =(lambdamax(band)-lambdamin(band))            ! in microns
    105122c        print*,band,lambdamin(band),dlambda,lambdamax(band)
    106123
    107 c changement de convention (signe) pour ksi,
    108 c et prise en compte de la largeur de bande (en cm-1):
     124c sign convention for ksi,
     125c and taking into account the wavelength width (in microns):
    109126         do mat=1,nbmat
    110127         do i=0,nlve+1
    111128           do j=0,nlve+1
    112               ksive(i,j,band,mat) = -ksive(i,j,band,mat)*dlambda
     129              ksive(i,j,band,mat) = +ksive(i,j,band,mat)*dlambda    ! in µm
    113130           enddo
    114131         enddo
    115132         enddo
    116 c calcul des coeff al et bl pour luminance Planck
     133c computing coeff al and bl for Planck luminance
    117134         al(band) = 2.*RHPLA*RCLUM*RCLUM/(lambda(band))**5.
    118 c cette luminance doit etre en W/m²/sr/µm pour correspondre au calcul
    119 c des ksi. Ici, elle est en W/m²/sr/m donc il faut mettre un facteur 1.e-6
     135c in W/m²/m
     136c We need W/m²/µm :
    120137     .                * 1.e-6
    121138         bl(band) = RHPLA*RCLUM/(RKBOL*lambda(band))

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

    r1017 r1301  
    4242c output
    4343
    44       REAL   PCOOL(klev) ! LONGWAVE COOLING (K/VENUSDAY) within each layer
     44      REAL   PCOOL(klev)  ! LONGWAVE COOLING (K/s) within each layer
    4545      REAL   PTOPLW       ! LONGWAVE FLUX AT T.O.A. (net, + vers le haut)
    4646      REAL   PSOLLW       ! LONGWAVE FLUX AT SURFACE (net, + vers le haut)
    … …  
    5151C* LOCAL VARIABLES:
    5252C
    53       real    dureejour
    54       parameter (dureejour=10.087e6)
    55      
    5653      integer i,j,p
    5754      real   zlnet(klev+1)    ! net thermal flux (W/m**2)
    … …  
    127124c     print*,PCOOL
    128125
    129       do j=1,klev
    130         PCOOL(j) = PCOOL(j)*dureejour ! K/Venusday
    131       enddo
    132 
    133126      return
    134127      end

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

    r1055 r1301  
    6969      REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx)
    7070      REAL albe(ngridmx),radsol(ngridmx),sollw(ngridmx)
    71       real solsw(ngridmx),dlw(ngridmx)
     71      real solsw(ngridmx),fder(ngridmx)
     72      real sollwdown(ngridmx),dlw(ngridmx)
    7273      REAL zmea(ngridmx), zstd(ngridmx)
    7374      REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx)
    … …  
    8485      real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
    8586      real albeS(ip1jmp1),radsolS(ip1jmp1),sollwS(ip1jmp1)
    86       real solswS(ip1jmp1),dlwS(ip1jmp1)
     87      real solswS(ip1jmp1),fderS(ip1jmp1)
     88      real sollwdownS(ip1jmp1),dlwS(ip1jmp1)
    8789      real zmeaS(ip1jmp1),zstdS(ip1jmp1),zsigS(ip1jmp1)
    8890      real zgamS(ip1jmp1),ztheS(ip1jmp1),zpicS(ip1jmp1)
    … …  
    116118      real, dimension(:,:),   allocatable :: albeold,radsolold
    117119      real, dimension(:,:),   allocatable :: sollwold,solswold
    118       real, dimension(:,:),   allocatable :: dlwold
     120      real, dimension(:,:),   allocatable :: fderold
     121      real, dimension(:,:),   allocatable :: dlwold,sollwdownold
    119122      real, dimension(:,:),   allocatable :: zmeaold,zstdold,zsigold
    120123      real, dimension(:,:),   allocatable :: zgamold,ztheold,zpicold
    … …  
    135138      integer, dimension(4) :: start,counter
    136139      REAL phisinverse(iip1,jjp1)  ! geopotentiel au sol avant inversion
    137       logical topoflag,albedoflag,razvitu
     140      logical topoflag,albedoflag,razvitu,razvitv
    138141      real    albedo
    139142     
    … …  
    318321      allocate(albeold(imold+1,jmold+1),radsolold(imold+1,jmold+1))
    319322      allocate(sollwold(imold+1,jmold+1),solswold(imold+1,jmold+1))
    320       allocate(dlwold(imold+1,jmold+1))
     323      allocate(fderold(imold+1,jmold+1))
     324      allocate(sollwdownold(imold+1,jmold+1),dlwold(imold+1,jmold+1))
    321325      allocate(zmeaold(imold+1,jmold+1),zstdold(imold+1,jmold+1))
    322326      allocate(zsigold(imold+1,jmold+1),zgamold(imold+1,jmold+1))
    … …  
    616620      ENDIF
    617621c
     622      ierr = NF_INQ_VARID (nid, "fder", nvarid)
     623      IF (ierr .NE. NF_NOERR) THEN
     624         PRINT*, "lect_start_archive: Le champ <fder> est absent"
     625         CALL abort
     626      ENDIF
     627#ifdef NC_DOUBLE
     628      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,fderold)
     629#else
     630      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,fderold)
     631#endif
     632      IF (ierr .NE. NF_NOERR) THEN
     633         PRINT*, "lect_start_archive: Lecture echouee pour <fder>"
     634         CALL abort
     635      ENDIF
     636c
    618637      ierr = NF_INQ_VARID (nid, "dlw", nvarid)
    619638      IF (ierr .NE. NF_NOERR) THEN
    … …  
    628647      IF (ierr .NE. NF_NOERR) THEN
    629648         PRINT*, "lect_start_archive: Lecture echouee pour <dlw>"
     649         CALL abort
     650      ENDIF
     651c
     652      ierr = NF_INQ_VARID (nid, "sollwdown", nvarid)
     653      IF (ierr .NE. NF_NOERR) THEN
     654         PRINT*, "lect_start_archive: Le champ <sollwdown> est absent"
     655         CALL abort
     656      ENDIF
     657#ifdef NC_DOUBLE
     658      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,sollwdownold)
     659#else
     660      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,sollwdownold)
     661#endif
     662      IF (ierr .NE. NF_NOERR) THEN
     663         PRINT*, "lect_start_archive: Lecture echouee pour <sollwdown>"
    630664         CALL abort
    631665      ENDIF
    … …  
    844878       topoflag = . FALSE .
    845879       CALL getin('topoflag',topoflag)
     880
     881        print*,zmeaold(2,1:10)
    846882
    847883       IF ( topoflag ) then
    … …  
    952988      call gr_dyn_fi (1,iip1,jjp1,ngridmx,solswS,solsw)
    953989
     990      call interp_horiz (fderold,fderS,imold,jmold,iim,jjm,1,
     991     &                   rlonuold,rlatvold,rlonu,rlatv)
     992      call gr_dyn_fi (1,iip1,jjp1,ngridmx,fderS,fder)
     993
    954994      call interp_horiz (dlwold,dlwS,imold,jmold,iim,jjm,1,
    955995     &                   rlonuold,rlatvold,rlonu,rlatv)
    956996      call gr_dyn_fi (1,iip1,jjp1,ngridmx,dlwS,dlw)
     997
     998      call interp_horiz (sollwdownold,sollwdownS,imold,jmold,iim,jjm,1,
     999     &                   rlonuold,rlatvold,rlonu,rlatv)
     1000      call gr_dyn_fi (1,iip1,jjp1,ngridmx,sollwdownS,sollwdown)
    9571001
    9581002      print*,"Nouvelles var physiques OK"
    … …  
    10441088       razvitu = . FALSE .
    10451089       CALL getin('razvitu',razvitu)
     1090       razvitv = . FALSE .
     1091       CALL getin('razvitv',razvitv)
    10461092
    10471093c calcul des champ de vent; passage en vent covariant
    … …  
    10531099      call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
    10541100     &                   rlonuold,rlatvold,rlonu,rlatv)
    1055       write (*,*) 'us ', us (1,2,1)   ! INFO
     1101      write (*,*) 'us ', us (1,2,1)    ! INFO
    10561102
    10571103      call interp_vert
    … …  
    10781124      write (*,*) 'ucov ', ucov (1,2,1)  ! INFO
    10791125c     write(48,*) 'ucov',ucov
     1126! Reseting v=0
     1127      if (razvitv) then
     1128           vnat(:,:,:) = 0.
     1129      endif
     1130      write (*,*) 'vnat ', vnat (1,2,1)    ! INFO
    10801131      do l=1,llm
    10811132        do j = 1, jjm
    … …  
    11141165
    11151166      call writerestartphy('restartphy.nc',tab_cntrl_fi,ngridmx,llm,
    1116      .           rlat,rlon,tsurf,tsoil,albe,solsw, sollw,dlw,
    1117      .           radsol,
     1167     .           rlat,rlon,tsurf,tsoil,albe,solsw, sollw,fder,dlw,
     1168     .           sollwdown,radsol,
    11181169     .           zmea,zstd,zsig,zgam,zthe,zpic,zval,
    11191170     .           t_fi)

  • trunk/LMDZ.VENUS/libf/phyvenus/phyetat0.F90

    r973 r1301  
    231231      PRINT*,'Rayonnement net au sol radsol:', xmin, xmax
    232232
    233 ! Lecture de l'orographie sous-maille si ok_orodr:
    234 
    235       if(ok_orodr) then
    236      
     233! Lecture de l'orographie sous-maille:
     234
    237235      CALL get_field("ZMEA",zmea,found)
    238236      IF (.not.found) THEN
    239237         PRINT*, 'phyetat0: Le champ <ZMEA> est absent'
    240          CALL abort
     238         PRINT*, 'mis a zero'
     239         zmea=0.
    241240      ENDIF
    242241      xmin = 1.0E+20
    … …  
    251250      IF (.not.found) THEN
    252251         PRINT*, 'phyetat0: Le champ <ZSTD> est absent'
    253          CALL abort
     252         PRINT*, 'mis a zero'
     253         zstd=0.
    254254      ENDIF
    255255      xmin = 1.0E+20
    … …  
    264264      IF (.not.found) THEN
    265265         PRINT*, 'phyetat0: Le champ <ZSIG> est absent'
    266          CALL abort
     266         PRINT*, 'mis a zero'
     267         zsig=0.
    267268      ENDIF
    268269      xmin = 1.0E+20
    … …  
    277278      IF (.not.found) THEN
    278279         PRINT*, 'phyetat0: Le champ <ZGAM> est absent'
    279          CALL abort
     280         PRINT*, 'mis a zero'
     281         zgam=0.
    280282      ENDIF
    281283      xmin = 1.0E+20
    … …  
    290292      IF (.not.found) THEN
    291293         PRINT*, 'phyetat0: Le champ <ZTHE> est absent'
    292          CALL abort
     294         PRINT*, 'mis a zero'
     295         zthe=0.
    293296      ENDIF
    294297      xmin = 1.0E+20
    … …  
    303306      IF (.not.found) THEN
    304307         PRINT*, 'phyetat0: Le champ <ZPIC> est absent'
    305          CALL abort
     308         PRINT*, 'mis a zero'
     309         zpic=0.
    306310      ENDIF
    307311      xmin = 1.0E+20
    … …  
    316320      IF (.not.found) THEN
    317321         PRINT*, 'phyetat0: Le champ <ZVAL> est absent'
    318          CALL abort
     322         PRINT*, 'mis a zero'
     323         zval=0.
    319324      ENDIF
    320325      xmin = 1.0E+20
    … …  
    325330      ENDDO
    326331      PRINT*,'OROGRAPHIE SOUS-MAILLE zval:', xmin, xmax
    327 
    328       else
    329          zmea = 0.
    330          zstd = 0.
    331          zsig = 0.
    332          zgam = 0.
    333          zthe = 0.
    334          zpic = 0.
    335          zval = 0.
    336 
    337       endif   ! fin test sur ok_orodr
    338332
    339333! Lecture de TANCIEN:

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

    r1160 r1301  
    77     .            paprs,pplay,ppk,pphi,pphis,presnivs,
    88     .            u,v,t,qx,
    9      .            omega,
     9     .            flxmw,
    1010     .            d_u, d_v, d_t, d_qx, d_ps)
    1111
    … …  
    4646c qx------input-R-mass mixing ratio traceurs (kg/kg)
    4747c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
    48 c omega---input-R-vitesse verticale en Pa/s
     48c flxmw---input-R-flux de masse vertical en kg/s
    4949c
    5050c d_u-----output-R-tendance physique de "u" (m/s/s)
    … …  
    125125      REAL d_t_dyn(klon,klev)
    126126
    127       REAL omega(klon,klev)
     127      REAL flxmw(klon,klev)
    128128
    129129      REAL d_u(klon,klev)
    … …  
    146146c Variables propres a la physique
    147147c
    148       REAL,save,allocatable :: rlev(:,:) ! altitude a chaque niveau (interface inferieure de la couche)
    149148      INTEGER,save :: itap        ! compteur pour la physique
    150149      REAL delp(klon,klev)        ! epaisseur d'une couche
     150      REAL omega(klon,klev)       ! vitesse verticale en Pa/s
     151
    151152     
    152153      INTEGER igwd,idx(klon),itest(klon)
    … …  
    242243c
    243244      REAL zphi(klon,klev)
    244 c
     245      REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
     246
    245247c Variables du changement
    246248c
    … …  
    362364c========================
    363365      IF (debut) THEN
    364          allocate(rlev(klon,klevp1))
    365366         allocate(source(klon,nqmax))
    366367
    … …  
    621622      ENDIF
    622623c====================================================================
     624
     625c Calcule de vitesse verticale a partir de flux de masse verticale
     626      DO k = 1, klev
     627       DO i = 1, klon
     628        omega(i,k) = RG*flxmw(i,k) / airephy(i)
     629       END DO
     630      END DO
     631
    623632c
    624633c Ajouter le geopotentiel du sol:
    … …  
    629638      ENDDO
    630639      ENDDO
     640
     641c   calcul du geopotentiel aux niveaux intercouches
     642c   ponderation des altitudes au niveau des couches en dp/p
     643
     644      DO k=1,klev
     645         DO i=1,klon
     646            zzlay(i,k)=zphi(i,k)/RG
     647         ENDDO
     648      ENDDO
     649      DO i=1,klon
     650         zzlev(i,1)=pphis(i)/RG
     651      ENDDO
     652      DO k=2,klev
     653         DO i=1,klon
     654            z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
     655            z2=(paprs(i,k)  +pplay(i,k))/(paprs(i,k)  -pplay(i,k))
     656            zzlev(i,k)=(z1*zzlay(i,k-1)+z2*zzlay(i,k))/(z1+z2)
     657         ENDDO
     658      ENDDO
     659      DO i=1,klon
     660         zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
     661      ENDDO
     662
    631663c====================================================================
    632664c
    … …  
    729761      fder = dlw
    730762
    731 c     print*,"radsol avant clmain=",radsol(klon/2)
    732 c     print*,"solsw avant clmain=",solsw(klon/2)
    733 c     print*,"sollw avant clmain=",sollw(klon/2)
    734 
    735763! ADAPTATION GCM POUR CP(T)
    736764
    … …  
    749777     s            ycoefh,yu1,yv1)
    750778
    751 c     print*,"radsol apres clmain=",radsol(klon/2)
    752 c     print*,"solsw apres clmain=",solsw(klon/2)
    753 c     print*,"sollw apres clmain=",sollw(klon/2)
    754 
    755779CXXX Incrementation des flux
    756780      DO i = 1, klon
    … …  
    770794      ENDDO
    771795      ENDDO
    772 c
    773 c        print*,"d_t_vdf1=",d_t_vdf(1,:)*dtime
    774 c        print*,"d_t_vdf2=",d_t_vdf(klon/2,:)*dtime
    775 c        print*,"d_t_vdf3=",d_t_vdf(klon,:)*dtime
    776 c        print*,"d_u_vdf=",d_u_vdf(klon/2,:)*dtime
    777 c        print*,"d_v_vdf=",d_v_vdf(klon/2,:)*dtime
    778796
    779797C TRACEURS
    … …  
    810828c Incrementer la temperature du sol
    811829c
    812 c     print*,'Tsol avant clmain:',ftsol(1)
    813830      DO i = 1, klon
    814831         ftsol(i) = ftsol(i) + d_ts(i)
    815832      ENDDO
    816 c     print*,'Tsol apres clmain:',ftsol(1)
    817833
    818834c Calculer la derive du flux infrarouge
    … …  
    878894         v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
    879895         d_v_ajs(:,:)= d_v_ajs(:,:)/dtime          ! (m/s)/s
    880       if (iflag_trac.eq.1) then
     896
     897         if (iflag_trac.eq.1) then
    881898           tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
    882899           d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime  ! /s
    883       endif
    884 
    885 c        print*,"d_t_ajs1=",d_t_ajs(1,:)*dtime
    886 c        print*,"d_t_ajs2=",d_t_ajs(klon/2,:)*dtime
    887 c        print*,"d_t_ajs3=",d_t_ajs(klon,:)*dtime
    888 c        print*,"d_u_ajs=",d_u_ajs(klon/2,:)*dtime
    889 c        print*,"d_v_ajs=",d_v_ajs(klon/2,:)*dtime
     900         endif
    890901
    891902      endif
    … …  
    919930c     PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
    920931           
    921 c     print*,"radsol avant radlwsw=",radsol(klon/2)
    922 c     print*,"solsw avant radlwsw=",solsw(klon/2)
    923 c     print*,"sollw avant radlwsw=",sollw(klon/2)
    924 c     print*,"avant radlwsw"
    925 
    926932      CALL radlwsw
    927      e            (dist, rmu0, fract,
     933     e            (dist, rmu0, fract, zzlev,
    928934     e             paprs, pplay,ftsol, t_seri,
    929935     s             heat,cool,radsol,
    … …  
    932938     s             lwnet, swnet)
    933939
    934 c     print*,"apres radlwsw"
    935 c     print*,"radsol apres radlwsw=",radsol(klon/2)
    936 c     print*,"solsw apres radlwsw=",solsw(klon/2)
    937 c     print*,"sollw apres radlwsw=",sollw(klon/2)
    938940      itaprad = 0
    939941      DO k = 1, klev
    940942      DO i = 1, klon
    941          dtrad(i,k) = heat(i,k)-cool(i,k)
    942          dtrad(i,k) = dtrad(i,k)/RDAY  !K/s
    943       ENDDO
    944       ENDDO
    945 c        print*,"dtrad1=",dtrad(1,:)*dtime
    946 c        print*,"dtrad2=",dtrad(klon/2,:)*dtime
    947 c        print*,"dtrad3=",dtrad(klon,:)*dtime
    948      
     943         dtrad(i,k) = heat(i,k)-cool(i,k)  !K/s
     944      ENDDO
     945      ENDDO
     946
    949947      ENDIF
    950948      itaprad = itaprad + 1

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

    r1160 r1301  
    7777      integer,parameter :: nbsrc=2,nblat=5,nblon=4
    7878      integer,parameter :: Nemiss=1   ! duree emission (Ed)
    79       real,save :: source_volcan(nbsrc)
     79      integer,save :: Nemiss(nbsrc)      ! duration emission (Ed)
     80      real,save :: source_volcan(nbsrc)  ! flux emission (kg/s)
    8081      real,save :: lat_volcan(nblat),lon_volcan(nblon)
    8182      real,save :: area_emiss(nblat,nblon)
    … …  
    131132         source_volcan(1) = 1.
    132133         source_volcan(2) = 1000.
     134c duration in Ed
     135         Nemiss(1) = 1
     136         Nemiss(2) = 10
    133137c localisation volcan
    134138         lat_volcan(1) =  70.
    … …  
    173177           do ilat  = 1,nblat
    174178            do ilon  = 1,nblon
    175              it=min(iemiss*ilat*ilon,nqtot)
     179             it=(iemiss-1)*nblat*nblon+(ilat-1)*nblon+ilon
     180             it=min(it,nqtot)
     181             deltatr(i,1,it) = 0.
    176182
    177183             if (i .eq. ig_volcan(ilat,ilon)) then
    178184
    179185c source appliquee pendant Nemiss Ed
    180                if (timesimu .lt. 86400.*Nemiss) then
     186               if (timesimu .lt. 86400.*Nemiss(iemiss)) then
    181187
    182188c source en kg/kg/s
    … …  
    185191           tr_seri(i,1,it) = tr_seri(i,1,it) + deltatr(i,1,it)*pdtphys
    186192
    187                else
    188            deltatr(i,1,it) = 0.
    189193               end if  ! duree emission
    190 
    191194             end if ! i localisation
    192195            end do

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

    r973 r1301  
    22! $Header: /home/cvsroot/LMDZ4/libf/phylmd/radlwsw.F,v 1.2 2004/10/27 10:14:46 lmdzadmin Exp $
    33!
    4       SUBROUTINE radlwsw(dist, rmu0, fract,
     4      SUBROUTINE radlwsw(dist, rmu0, fract, zzlev,
    55     .                  paprs, pplay,tsol, t,
    66     .                  heat,cool,radsol,
    … …  
    1717c fract----input-R- duree d'ensoleillement normalisee
    1818c solaire--input-R- constante solaire (W/m**2) (dans clesphys.h)
     19c zzlev----input-R- altitude a inter-couche (m)
    1920c paprs----input-R- pression a inter-couche (Pa)
    2021c pplay----input-R- pression au milieu de couche (Pa)
    2122c tsol-----input-R- temperature du sol (en K)
    2223c t--------input-R- temperature (K)
    23 c heat-----output-R- echauffement atmospherique (visible) (K/jour)
    24 c cool-----output-R- refroidissement dans l'IR (K/jour)
     24c heat-----output-R- echauffement atmospherique (visible) (K/s)
     25c cool-----output-R- refroidissement dans l'IR (K/s)
    2526c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
    2627c topsw----output-R- flux solaire net au sommet de l'atm. (+ vers le bas)
    … …  
    3637c   S. Lebonnois    20/12/2006
    3738c   corrections     13/07/2007
     39c   New ksi matrix: possibility of different cloud model fct of lat   05/2014
    3840
    3941c======================================================================
    … …  
    4648#include "clesphys.h"
    4749#include "comcstVE.h"
    48 c
     50
     51!===========
     52! Arguments
     53!===========
    4954      real rmu0(klon), fract(klon), dist
    50 c
     55
     56      REAL zzlev(klon,klev+1)
    5157      real paprs(klon,klev+1), pplay(klon,klev)
    5258      real tsol(klon)
    … …  
    5763      real sollwdown(klon)
    5864      REAL swnet(klon,klev+1),lwnet(klon,klev+1)
    59 c
     65
     66!===========
     67! Local
     68!===========
    6069      INTEGER k, kk, i, j, band
    61 c
     70
    6271      REAL   PPB(klev+1)
    63 c
     72
    6473      REAL   zfract, zrmu0
    65 c
     74
    6675      REAL   zheat(klev), zcool(klev)
    67       real  temp(klev)
     76      real   temp(klev),znivs(klev+1)
    6877      REAL   ZFSNET(klev+1),ZFLNET(klev+1)
    6978      REAL   ztopsw, ztoplw
    … …  
    7382cIM END
    7483      real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file
    75       real,save,allocatable :: ztop(:) ! in km
    7684
    7785      real    psi(0:klev+1,0:klev+1)
    7886      real    deltapsi(0:klev+1,0:klev+1)
    79       real    latdeg
    8087      real    pt0(0:klev+1)
    8188      real    bplck(0:klev+1,nnuve)    ! Planck luminances in table layers
    82       real    y(0:klev,nnuve)          ! intermediaire Planck
    83       real    zdblay(0:klev+1,nnuve)   ! gradient en temperature de planck
    84       integer mat,mat0
     89      real    y(0:klev,nnuve)          ! temporary variable for Planck
     90      real    zdblay(0:klev+1,nnuve)   ! temperature gradient of planck function
     91      integer mat0,lat,ips,isza,ips0,isza0
    8592      real    factp,factz,ksi
    8693
    … …  
    98105        allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat))
    99106        call load_ksi(ksive)
    100 c ---------- ztop --------------
    101         allocate(ztop(klon))
    102         DO i = 1, klon
    103              ztop(i) = 70.
    104         ENDDO !i
    105 c ztop: d'apres fit à figure 16 du papier Zavosa et al (tmp) traitant des
    106 c       donnees Venera
    107 c       DO i = 1, klon
    108 c         latdeg = abs(rlatd(i))
    109 c         if (latdeg.lt.15) then
    110 c            ztop(i) = 70.
    111 c         elseif (latdeg.lt.50) then
    112 c            ztop(i) = 63.95+6*cos((latdeg-15)*RPI/2./50.)
    113 c         else
    114 c            ztop(i) = min(63.95+6*cos((latdeg-15)*RPI/2./50.),
    115 c    .                     63.95-5.9*sin((latdeg-60)*RPI/2/30))
    116 c         endif
    117 c       print*,'lat(',i,')=',latdeg,'  ztop=',ztop(i)
    118 c       ENDDO !i
    119 c ---------- ztop --------------
    120107
    121108      endif ! firstcall
    … …  
    191178c --------------------------
    192179       mat0 = 0
    193        do mat=1,nbmat-nbztopve
    194          if (  (psurfve(mat).ge.paprs(j,1))
    195      .    .and.(psurfve(mat+nbztopve).lt.paprs(j,1))
    196      .    .and.(ztopve(mat).lt.ztop(j))
    197      .    .and.(ztopve(mat+1).ge.ztop(j)) ) then
    198               mat0  = mat
    199 c             print*,'ig=',j,'  mat0=',mat
    200               factp = (paprs(j,1)           -psurfve(mat))
    201      .               /(psurfve(mat+nbztopve)-psurfve(mat))
    202               factz = (ztop(j)      -ztopve(mat))
    203      .               /(ztopve(mat+1)-ztopve(mat))
     180       if (nlatve.eq.1) then   ! clouds are taken as uniform
     181         lat=1
     182       else
     183         if (abs(rlatd(j)).le.50.) then
     184           lat=1
     185         elseif (abs(rlatd(j)).le.60.) then
     186           lat=2
     187         elseif (abs(rlatd(j)).le.70.) then
     188           lat=3
     189         elseif (abs(rlatd(j)).le.80.) then
     190           lat=4
     191         else
     192           lat=5
     193         endif
     194       endif
     195       
     196       ips0=0
     197       do ips=1,nbpsve(lat)-1
     198         if (  (psurfve(ips,lat).ge.paprs(j,1))
     199     .    .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then
     200              ips0  = ips
     201c             print*,'ig=',j,'  ips0=',ips
     202              factp = (paprs(j,1)         -psurfve(ips0,lat))
     203     .               /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
    204204              exit
    205205         endif
    206206       enddo
    207        if (mat0.eq.0) then
     207       isza0=0
     208       if (nbszave(lat).gt.1) then
     209        do isza=1,nbszave(lat)-1
     210         if (  (szave(isza,lat).ge.zrmu0)
     211     .    .and.(szave(isza+1,lat).lt.zrmu0) ) then
     212              isza0  = isza
     213c             print*,'ig=',j,'  isza0=',isza
     214              factz = (zrmu0             -szave(isza0,lat))
     215     .               /(szave(isza0+1,lat)-szave(isza0,lat))
     216              exit
     217         endif
     218        enddo
     219       else            ! Only one sza, no interpolation
     220        isza0=-99
     221       endif
     222       
     223       if ((ips0.eq.0).or.(isza0.eq.0)) then
    208224         write(*,*) 'Finding the right matrix in radlwsw'
    209          print*,'Probleme pour interpolation au point ig=',j
    210          print*,'psurf = ',paprs(j,1),' ztop = ',ztop(j)
     225         print*,'Interpolation problem, grid point ig=',j
     226         print*,'psurf = ',paprs(j,1),' mu0 = ',zrmu0
    211227         stop
    212228       endif
    213 
     229       
     230       if (isza0.eq.-99) then
     231         mat0 = indexve(lat)+ips0
     232       else
     233         mat0 = indexve(lat)+(isza0-1)*nbpsve(lat)+ips0
     234       endif
     235       
    214236c interpolation of ksi and computation of psi,deltapsi
    215237c ----------------------------------------------------
    … …  
    217239        do k=0,klev+1
    218240         do i=0,klev+1
    219           ksi = ksive(i,k,band,mat0)*(1-factz)*(1-factp)
    220      .         +ksive(i,k,band,mat0+1)*factz  *(1-factp)
    221      .         +ksive(i,k,band,mat0+nbztopve)*(1-factz)*factp
    222      .         +ksive(i,k,band,mat0+nbztopve+1)*factz  *factp
     241          if (isza0.eq.-99) then
     242           ksi = ksive(i,k,band,mat0)*(1-factp)
     243     .          +ksive(i,k,band,mat0+1)*factp
     244          else
     245           ksi = ksive(i,k,band,mat0)*(1-factp)*(1-factz)
     246     .          +ksive(i,k,band,mat0+1)*factp  *(1-factz)
     247     .          +ksive(i,k,band,mat0+nbpsve(lat))*(1-factp)*factz
     248     .          +ksive(i,k,band,mat0+nbpsve(lat)+1)*factp  *factz
     249          endif
     250     
    223251          psi(i,k) = psi(i,k) +
    224      .               ksi*(bplck(i,band)-bplck(k,band))
    225           deltapsi(i,k) = deltapsi(i,k) + ksi*zdblay(i,band)
     252     .               RPI*ksi*(bplck(i,band)-bplck(k,band))
     253          deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
    226254         enddo
    227255        enddo
    … …  
    241269c SW call
    242270c---------
     271      znivs=zzlev(j,:)
     272c      CALL SW_venus_ve(zrmu0, zfract,
     273c     S        PPB,temp,znivs,
     274c     S        zheat,
     275c     S        ztopsw,zsolsw,ZFSNET)
     276
    243277      CALL SW_venus_dc(zrmu0, zfract,
    244278     S        PPB,temp,
    … …  
    283317c     j = 1
    284318c     print*,'mu0=',rmu0(j)
    285 c     print*,'   net flux vis   HEAT(K/day)'
     319c     print*,'   net flux vis   HEAT(K/Eday)'
    286320c     do k=1,klev
    287 c     print*,k,ZFSNET(k),heat(j,k)*8.56548e-3
     321c     print*,k,ZFSNET(k),heat(j,k)*86400.
    288322c     enddo
    289 c     print*,'   net flux IR    COOL(K/day)'
     323c     print*,'   net flux IR    COOL(K/Eday)'
    290324c     do k=1,klev
    291 c     print*,k,ZFLNET(k),cool(j,k)*8.56548e-3
     325c     print*,k,ZFLNET(k),cool(j,k)*86400.
    292326c     enddo
    293327

  • trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.NewtonCool

    r892 r1301  
    22! $Header: /home/cvsroot/LMDZ4/libf/phylmd/radlwsw.F,v 1.2 2004/10/27 10:14:46 lmdzadmin Exp $
    33!
    4       SUBROUTINE radlwsw(dist, rmu0, fract,
    5      .                  paprs, pplay,tsol, t,
    6      .                  heat,cool,radsol,
    7      .                  topsw,toplw,solsw,sollw,
    8      .                  sollwdown,
    9      .                  lwnet, swnet)
    10 c     
     4      SUBROUTINE radlwsw(dist, rmu0, fract, zzlev,
     5     .                  paprs, pplay,tsol, pt, nq, nmicro, pq,qaer)
     6     
    117c======================================================================
    128c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
    … …  
    1612c rmu0-----input-R- cosinus de l'angle zenithal
    1713c fract----input-R- duree d'ensoleillement normalisee
    18 c solaire--input-R- constante solaire (W/m**2) (dans clesphys.h)
    1914c paprs----input-R- pression a inter-couche (Pa)
    2015c pplay----input-R- pression au milieu de couche (Pa)
    2116c tsol-----input-R- temperature du sol (en K)
    22 c t--------input-R- temperature (K)
    23 c heat-----output-R- echauffement atmospherique (visible) (K/jour)
    24 c cool-----output-R- refroidissement dans l'IR (K/jour)
    25 c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
    26 c topsw----output-R- flux solaire net au sommet de l'atm. (+ vers le bas)
    27 c toplw----output-R- ray. IR net au sommet de l'atmosphere (+ vers le haut)
    28 c solsw----output-R- flux solaire net a la surface (+ vers le bas)
    29 c sollw----output-R- ray. IR net a la surface (+ vers le bas)
    30 c sollwdown-output-R- ray. IR descendant a la surface (+ vers le bas)
    31 c lwnet____output-R- flux IR net (+ vers le haut)
    32 c swnet____output-R- flux solaire net (+ vers le bas)
     17c pt-------input-R- temperature (K)
    3318c
    3419     
    3520c   S. Lebonnois    12/04/2007
    3621c  VERSION NEWTONIAN COOLING pour Venus (no diurnal cycle)
     22c  update 01/2014
    3723
    3824c======================================================================
    3925      use dimphy
    4026      USE comgeomphy
     27      USE phys_state_var_mod, only: falbe,heat,cool,radsol,
     28     .      topsw,toplw,solsw,sollw,sollwdown,lwnet,swnet
     29      USE write_field_phy
    4130      IMPLICIT none
    4231#include "dimensions.h"
    4332#include "YOMCST.h"
    4433#include "clesphys.h"
    45 c
     34
     35c ARGUMENTS
     36      INTEGER nq,nmicro
    4637      real rmu0(klon), fract(klon), dist
    47 c
    48       real paprs(klon,klev+1), pplay(klon,klev)
     38 
     39      real zzlev(klon,klev+1),paprs(klon,klev+1), pplay(klon,klev)
    4940      real tsol(klon)
    50       real t(klon,klev)
    51       real heat(klon,klev), cool(klon,klev)
    52       real radsol(klon), topsw(klon), toplw(klon)
    53       real solsw(klon), sollw(klon)
    54       real sollwdown(klon)
    55       REAL swnet(klon,klev+1),lwnet(klon,klev+1)
    56 c
     41      real pt(klon,klev)
     42      real pq(klon,klev,nq)
     43      REAL qaer(klon,klev,nq)
     44 
     45c LOCAL VARIABLES
    5746      INTEGER i,j,k
    5847      integer    nlevCLee,level
    … …  
    6453      real   ztemp,zdt,fact
    6554      real   dTsdt(klev),zt_eq(klon,klev)
    66       real   dureejour
    67       parameter (dureejour=10.087e6)
    6855      save   zt_eq
    6956     
    … …  
    144131
    145132        DO k = 1, klev
    146          heat (j,k) = dTsdt(k)*dureejour ! K/Venusday
     133         heat (j,k) = dTsdt(k)    ! K/s
    147134         cool (j,k) = 0.
    148135        ENDDO

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

    r1300 r1301  
    66      USE phys_state_var_mod
    77      use cpdet_mod, only: ini_cpdet
     8      use moyzon_mod, only: tmoy
     9
    810      IMPLICIT NONE
    911
    … …  
    288290        print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
    289291      ENDDO
     292
     293      allocate(tmoy(llm))
     294      tmoy(:)=temp(:)
    290295     
    291296c     temperature du sous-sol
    … …  
    465470#include "../dyn3d_common/disvert_noterre.F"
    466471#include "../dyn3d/abort_gcm.F"
    467 !#include "../dyn3d/dump2d.F"
    468 
     472

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

    r808 r1301  
    77     .            rlat,rlon, tsol,tsoil,
    88     .           albe, solsw, sollw,
    9      .           fder,radsol,
     9     .           fder,dlw, sollwdown, radsol,
    1010     .    zmea, zstd, zsig, zgam, zthe, zpic, zval,
    1111     .           tabcntr0)
    … …  
    2929      REAL tsoil(ngridmx,nsoilmx)
    3030      REAL albe(ngridmx)
    31 cIM BEG alblw
    32       REAL alblw(ngridmx)
    33 cIM END alblw
    3431      REAL radsol(ngridmx)
    3532      REAL sollw(ngridmx)
    3633      real solsw(ngridmx)
    3734      real fder(ngridmx)
     35      real dlw(ngridmx)
     36      real sollwdown(ngridmx)
    3837      REAL zmea(ngridmx), zstd(ngridmx)
    3938      REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx)
    … …  
    275274
    276275c
     276c Lecture derive des flux IR:
     277c
     278      ierr = NF_INQ_VARID (nid, "dlw", nvarid)
     279      IF (ierr.NE.NF_NOERR) THEN
     280         PRINT*, 'phyetat0: Le champ <dlw> est absent'
     281         PRINT*, 'mis a zero'
     282         dlw = 0.
     283      ELSE
     284#ifdef NC_DOUBLE
     285        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, dlw)
     286#else
     287        ierr = NF_GET_VAR_REAL(nid, nvarid, dlw)
     288#endif
     289        IF (ierr.NE.NF_NOERR) THEN
     290          PRINT*, 'phyetat0: Lecture echouee pour <dlw>'
     291          CALL abort
     292        ENDIF
     293      ENDIF
     294      xmin = 1.0E+20
     295      xmax = -1.0E+20
     296      DO i = 1, ngridmx
     297         xmin = MIN(dlw(i),xmin)
     298         xmax = MAX(dlw(i),xmax)
     299      ENDDO
     300      PRINT*,'Derive des flux IR dlw:', xmin, xmax
     301
     302c
     303c Lecture rayonnement IR vers le bas au sol:
     304c
     305      ierr = NF_INQ_VARID (nid, "sollwdown", nvarid)
     306      IF (ierr.NE.NF_NOERR) THEN
     307         PRINT*, 'phyetat0: Le champ <sollwdown> est absent'
     308         PRINT*, 'mis a zero'
     309         sollwdown = 0.
     310      ELSE
     311#ifdef NC_DOUBLE
     312        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, sollwdown)
     313#else
     314        ierr = NF_GET_VAR_REAL(nid, nvarid, sollwdown)
     315#endif
     316        IF (ierr.NE.NF_NOERR) THEN
     317          PRINT*, 'phyetat0: Lecture echouee pour <sollwdown>'
     318          CALL abort
     319        ENDIF
     320      ENDIF
     321      xmin = 1.0E+20
     322      xmax = -1.0E+20
     323      DO i = 1, ngridmx
     324         xmin = MIN(sollwdown(i),xmin)
     325         xmax = MAX(sollwdown(i),xmax)
     326      ENDDO
     327      PRINT*,'Flux IR vers le bas au sol sollwdown:', xmin, xmax
     328
     329c
    277330c Lecture du rayonnement net au sol:
    278331c

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

    r1055 r1301  
    1818      USE infotrac
    1919      USE control_mod
    20       use cpdet_mod, only: tpot2t
     20      use cpdet_mod, only: tpot2t,ini_cpdet
    2121
    2222      implicit none
    … …  
    4242c -----------------------------
    4343      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    44       REAL teta(ip1jmp1,llm)                    ! temperature potentielle
     44      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
    4545      REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
    4646      REAL pks(ip1jmp1)                      ! exner (f pour filtre)
    … …  
    5151      REAL masse(ip1jmp1,llm)                ! masse de l'atmosphere
    5252      REAL ps(ip1jmp1)                       ! pression au sol
    53       REAL p3d(iip1, jjp1, llm+1)            ! pression aux interfaces
     53      REAL p3d(iip1,jjp1,llm+1)              ! pression aux interfaces
    5454     
    5555c Variable Physiques (grille physique)
    … …  
    6262      REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx)
    6363      REAL albe(ngridmx),radsol(ngridmx),sollw(ngridmx)
    64       real solsw(ngridmx),dlw(ngridmx)
     64      real solsw(ngridmx),fder(ngridmx)
     65      real sollwdown(ngridmx),dlw(ngridmx)
    6566      REAL zmea(ngridmx), zstd(ngridmx)
    6667      REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx)
    … …  
    6970      INTEGER start,length
    7071      PARAMETER (length = 100)
    71       REAL tab_cntrl_fi(length) ! tableau des parametres de startfi
     72      REAL tab_cntrl_fi(length)  ! tableau des parametres de startfi
    7273      REAL tab_cntrl_dyn(length) ! tableau des parametres de start
    7374      INTEGER*4 day_ini_fi
    … …  
    7980      real rlatS(ip1jmp1),rlonS(ip1jmp1)
    8081      real albeS(ip1jmp1),radsolS(ip1jmp1),sollwS(ip1jmp1)
    81       real solswS(ip1jmp1),dlwS(ip1jmp1)
     82      real solswS(ip1jmp1),fderS(ip1jmp1)
     83      real dlwS(ip1jmp1),sollwdownS(ip1jmp1)
    8284      real zmeaS(ip1jmp1),zstdS(ip1jmp1),zsigS(ip1jmp1)
    8385      real zgamS(ip1jmp1),ztheS(ip1jmp1),zpicS(ip1jmp1)
    … …  
    132134     .       rlat,rlon,tsurf,tsoil,
    133135     .       albe, solsw, sollw,
    134      .       dlw,radsol,
     136     .       fder,dlw,sollwdown,radsol,
    135137     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,
    136138     .       tab_cntrl_fi)
    … …  
    140142c-----------------------------------------------------------------------
    141143
     144      CALL conf_gcm( 99, .TRUE. )
    142145      call iniconst
    143146      call inigeom
    144147      call inifilr
     148      call ini_cpdet
     149
    145150      CALL pression(ip1jmp1, ap, bp, ps, p3d)
    146151         if (disvert_type==1) then
    … …  
    198203      call gr_fi_dyn(1,ngridmx,iip1,jjp1,sollw,sollwS)
    199204      call gr_fi_dyn(1,ngridmx,iip1,jjp1,solsw,solswS)
     205      call gr_fi_dyn(1,ngridmx,iip1,jjp1,fder,fderS)
    200206      call gr_fi_dyn(1,ngridmx,iip1,jjp1,dlw,dlwS)
     207      call gr_fi_dyn(1,ngridmx,iip1,jjp1,sollwdown,sollwdownS)
    201208      call gr_fi_dyn(1,ngridmx,iip1,jjp1,zmea,zmeaS)
    202209      call gr_fi_dyn(1,ngridmx,iip1,jjp1,zstd,zstdS)
    … …  
    295302      call write_archive(nid,'solsw',
    296303     .             'SW flux at surface','W m-2',2,solswS)
     304      call write_archive(nid,'fder','derive','?',2,fderS)
    297305      call write_archive(nid,'dlw','LW derive','?',2,dlwS)
     306      call write_archive(nid,'sollwdown',
     307     .             'LW dwn flux at surface','?',2,sollwdownS)
    298308      call write_archive(nid,'zmea','param oro sous-maille','m',2,zmeaS)
    299309      call write_archive(nid,'zstd','param oro sous-maille','m',2,zstdS)

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

    r1017 r1301  
    4040c output
    4141
    42       REAL   PHEAT(klev) ! SHORTWAVE HEATING (K/VENUSDAY) within each layer
     42      REAL   PHEAT(klev)  ! SHORTWAVE HEATING (K/s) within each layer
    4343      REAL   PTOPSW       ! SHORTWAVE FLUX AT T.O.A. (net)
    4444      REAL   PSOLSW       ! SHORTWAVE FLUX AT SURFACE (net)
    … …  
    4949C
    5050      integer nldc,nszadc
    51       real    dureejour
    5251      parameter (nldc=49)  ! fichiers Crisp
    5352      parameter (nszadc=8) ! fichiers Crisp
    54       parameter (dureejour=10.087e6)
    5553     
    5654      integer i,j,nsza,nsza0,nl0
    … …  
    234232        PHEAT(j) = (ZFSNET(j+1)-ZFSNET(j))
    235233     .            *RG/cpdet(pt(j)) / ((PPB(j)-PPB(j+1))*1.e5)
    236         PHEAT(j) = PHEAT(j)*dureejour  ! K/venus_day
    237234      enddo
    238235

  • trunk/LMDZ.VENUS/libf/phyvenus/write_histday.h

    r902 r1301  
    8686c en K/s     
    8787      call histwrite_phy(nid_day,.false.,"dtajs",itau_w,d_t_ajs)
    88 c K/day ==> K/s
    89       call histwrite_phy(nid_day,.false.,"dtswr",itau_w,heat/RDAY)
    90 c K/day ==> K/s     
    91       call histwrite_phy(nid_day,.false.,"dtlwr",itau_w,-1.*cool/RDAY)
     88c en K/s     
     89      call histwrite_phy(nid_day,.false.,"dtswr",itau_w,heat)
     90c en K/s     
     91      call histwrite_phy(nid_day,.false.,"dtlwr",itau_w,-1.*cool)
    9292c en K/s     
    9393c     call histwrite_phy(nid_day,.false.,"dtec",itau_w,d_t_ec)

  • trunk/LMDZ.VENUS/libf/phyvenus/write_histins.h

    r902 r1301  
    8686c en K/s     
    8787      call histwrite_phy(nid_ins,.false.,"dtajs",itau_w,d_t_ajs)
    88 c K/day ==> K/s
    89       call histwrite_phy(nid_ins,.false.,"dtswr",itau_w,heat/RDAY)
    90 c K/day ==> K/s     
    91       call histwrite_phy(nid_ins,.false.,"dtlwr",itau_w,-1.*cool/RDAY)
     88c en K/s     
     89      call histwrite_phy(nid_ins,.false.,"dtswr",itau_w,heat)
     90c en K/s     
     91      call histwrite_phy(nid_ins,.false.,"dtlwr",itau_w,-1.*cool)
    9292c en K/s     
    9393c     call histwrite_phy(nid_ins,.false.,"dtec",itau_w,d_t_ec)

  • trunk/LMDZ.VENUS/libf/phyvenus/write_histmth.h

    r902 r1301  
    8686c en K/s     
    8787      call histwrite_phy(nid_mth,.false.,"dtajs",itau_w,d_t_ajs)
    88 c K/day ==> K/s
    89       call histwrite_phy(nid_mth,.false.,"dtswr",itau_w,heat/RDAY)
    90 c K/day ==> K/s     
    91       call histwrite_phy(nid_mth,.false.,"dtlwr",itau_w,-1.*cool/RDAY)
     88c en K/s     
     89      call histwrite_phy(nid_mth,.false.,"dtswr",itau_w,heat)
     90c en K/s     
     91      call histwrite_phy(nid_mth,.false.,"dtlwr",itau_w,-1.*cool)
    9292c en K/s     
    9393c     call histwrite_phy(nid_mth,.false.,"dtec",itau_w,d_t_ec)

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

    r808 r1301  
    22     .           rlat,rlon, tsol,tsoil,
    33     .           albedo,
    4      .           solsw, sollw,fder,
    5      .           radsol,
     4     .           solsw, sollw,fder,dlw,
     5     .           sollwdown,radsol,
    66     .    zmea, zstd, zsig, zgam, zthe, zpic, zval,
    77     .           t_ancien)
    … …  
    2727      real sollw(klon)
    2828      real fder(klon)
     29      real dlw(klon)
     30      real sollwdown(klon)
    2931      REAL radsol(klon)
    3032      REAL zmea(klon), zstd(klon)
    … …  
    182184#endif
    183185      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 27,
    184      .                        "Rayonnement IF a la surface")
     186     .                        "Rayonnement IR a la surface")
    185187      ierr = NF_ENDDEF(nid)
    186188#ifdef NC_DOUBLE
    … …  
    207209      ierr = NF_REDEF (nid)
    208210#ifdef NC_DOUBLE
     211      ierr = NF_DEF_VAR (nid, "dlw", NF_DOUBLE, 1, idim2,nvarid)
     212#else
     213      ierr = NF_DEF_VAR (nid, "dlw", NF_FLOAT, 1, idim2,nvarid)
     214#endif
     215      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 14,
     216     .                        "Derivee flux IR")
     217      ierr = NF_ENDDEF(nid)
     218#ifdef NC_DOUBLE
     219      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,dlw)
     220#else
     221      ierr = NF_PUT_VAR_REAL (nid,nvarid,dlw)
     222#endif
     223c
     224      ierr = NF_REDEF (nid)
     225#ifdef NC_DOUBLE
     226      ierr = NF_DEF_VAR (nid, "sollwdown", NF_DOUBLE, 1, idim2,nvarid)
     227#else
     228      ierr = NF_DEF_VAR (nid, "sollwdown", NF_FLOAT, 1, idim2,nvarid)
     229#endif
     230      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 14,
     231     .                        "Flux IR vers le bas a la surface")
     232      ierr = NF_ENDDEF(nid)
     233#ifdef NC_DOUBLE
     234      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,sollwdown)
     235#else
     236      ierr = NF_PUT_VAR_REAL (nid,nvarid,sollwdown)
     237#endif
     238c
     239      ierr = NF_REDEF (nid)
     240#ifdef NC_DOUBLE
    209241      ierr = NF_DEF_VAR (nid, "RADS", NF_DOUBLE, 1, idim2,nvarid)
    210242#else

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

    r892 r1301  
    2121c lat------INPUT : latitude en degres
    2222c long-----INPUT : longitude en degres
    23 c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
    24 c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
     23c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad/RDAY
     24c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad/RDAY
    2525c================================================================
    2626      use dimphy
    … …  
    9898      IF (omega1.LE.omega2) THEN  !--on est dans la meme journee locale
    9999c
    100       IF (omega2.LE.-omega .OR. omega1.GE.omega
    101      .                     .OR. omega.LT.1e-5) THEN   !--nuit
    102          frac(i)=0.0
    103          pmu0(i)=0.0
    104       ELSE                                              !--jour+nuit/jour
     100      IF (omega2.LE.-omega .OR. omega1.GE.omega      !--nuit
     101     .           .OR. omega.LT.1e-5) THEN            !--nuit polaire
     102        frac(i)=0.0
     103        pmu0(i)=SIN(latr)*SIN(lat_sun) +
     104     .          COS(latr)*COS(lat_sun)*
     105     .          (SIN(omega2)-SIN(omega1))/
     106     .          (omega2-omega1)       
     107      ELSE                                           !--jour+nuit/jour
    105108        omegadeb=MAX(-omega,omega1)
    106109        omegafin=MIN(omega,omega2)
    … …  
    117120      IF (omega1.GE.omega) THEN  !--nuit
    118121         zfrac1=0.0
    119          z1_mu =0.0
     122         z1_mu =SIN(latr)*SIN(lat_sun) +
     123     .          COS(latr)*COS(lat_sun)*
     124     .          (-SIN(omega1))/       
     125     .          (RPI-omega1)
    120126      ELSE                       !--jour+nuit
    121127        omegadeb=MAX(-omega,omega1)
    … …  
    130136      IF (omega2.LE.-omega) THEN   !--nuit
    131137         zfrac2=0.0
    132          z2_mu =0.0
     138         z2_mu =SIN(latr)*SIN(lat_sun) +
     139     .          COS(latr)*COS(lat_sun)*
     140     .          (SIN(omega2))/       
     141     .          (omega2+RPI)
    133142      ELSE                         !--jour+nuit
    134143         omegadeb=-omega
    … …  
    143152c-----------------------moyenne
    144153      frac(i)=(zfrac1+zfrac2)/(omega2+2*RPI-omega1)
    145       pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10)
     154      if (frac(i).ne.0.) then
     155       pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10)
     156      else
     157       pmu0(i)=((RPI-omega1)*z1_mu+(omega2+RPI)*z2_mu)/
     158     .                   (omega2+2*RPI-omega1)
     159      endif
    146160c
    147161      ENDIF   !---comparaison omega1 et omega2
Note: See TracChangeset for help on using the changeset viewer.