Changeset 652


Ignore:
Timestamp:
May 7, 2012, 4:01:14 PM (13 years ago)
Author:
acolaitis
Message:

Thermals:

  • minor changes

Diffusion with Yamada4:

  • corrected numerical oscillation in daytime first-level temperature and momentum predictions through the implementation of a sub-timestep (40 sub-steps) of one loop. (cheap in cpu time). This loop is the one that computes the tke at t+dt from the tke at t.
  • corrected numerical oscillations in nighttime first-level temperature and momentum predictions through the use of a minimum mixing coefficient, computed as in Holtslag and Boville (93), tuned for Mars.
  • yamada4 now takes into account the diffusion of tke (iflag_pbl=9), as it tends to smooth numerical oscillations. Results are very similar to those obtained without diffusion.
Location:
trunk/LMDZ.MARS/libf/phymars
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90

    r628 r652  
    33!
    44      SUBROUTINE calltherm_interface (firstcall, &
    5      & long,lati,zzlev,zzlay, &
     5     & zzlev,zzlay, &
    66     & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, &
    77     & pplay,pplev,pphi,zpopsk, &
    … …  
    2121!--------------------------------------------------------
    2222
     23!      REAL, INTENT(IN) :: long(ngridmx),lati(ngridmx)
    2324      REAL, INTENT(IN) :: ptimestep
    2425      REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1),pplay(ngridmx,nlayermx)
    … …  
    3233      REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx),pdt(ngridmx,nlayermx)
    3334      REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1)
    34       REAL, INTENT(IN) :: long(ngridmx),lati(ngridmx)
    3535      REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx)
    3636
    … …  
    363363      if (dtke_thermals) then
    364364      detrmod(:,:)=0.
    365       ndt=1
     365      ndt=10
    366366      do l=1,zlmax
    367367         do ig=1,ngridmx

  • trunk/LMDZ.MARS/libf/phymars/physiq.F

    r636 r652  
    846846 
    847847        call calltherm_interface(firstcall,
    848      $ long,lati,zzlev,zzlay,
     848     $ zzlev,zzlay,
    849849     $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
    850850     $ pplay,pplev,pphi,zpopsk,

  • trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90

    r648 r652  
    857857            if (l.le.lmax(ig)) then
    858858                if (zw2(ig,l).lt.0.)then
    859                   print*,'pb2 zw2<0'
     859!                  print*,'pb2 zw2<0',zw2(ig,l)
     860                  zw2(ig,l)=0.
    860861                endif
    861862                zw2(ig,l)=sqrt(zw2(ig,l))

  • trunk/LMDZ.MARS/libf/phymars/vdifc.F

    r626 r652  
    384384      CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
    385385     s   ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true,pq2,zkv,zkh,zkq,ust
    386      s   ,8)
     386     s   ,9)
    387387
    388388      ENDIF

  • trunk/LMDZ.MARS/libf/phymars/yamada4.F

    r554 r652  
    120120      REAL teta(ngrid,nlay)
    121121      REAL pq(ngrid,nlay,nqmx)
     122      REAL kminfact
     123      INTEGER i
     124      REAL ztimestep
     125      INTEGER, SAVE :: ndt
    122126
    123127      nlev=nlay+1
    … …  
    146150        endif
    147151      allocate(l0(ngrid))
     152      ndt=ceiling(3840./(3699.*24./dt))
    148153      firstcall=.false.
    149154      endif !of if firstcall
    … …  
    153158c.......................................................................
    154159
    155       if (ico2.ne.0) then
    156 !     Special case if one of the tracers is CO2 gas
    157          DO k=1,nlay
    158            DO ig=1,ngrid
    159             teta(ig,k) = phc(ig,k)*(A*pq(ig,k,ico2)+B)
    160            ENDDO
    161          ENDDO
    162        else
     160!      if (ico2.ne.0) then
     161!!     Special case if one of the tracers is CO2 gas
     162!         DO k=1,nlay
     163!           DO ig=1,ngrid
     164!            teta(ig,k) = phc(ig,k)*(A*pq(ig,k,ico2)+B)
     165!           ENDDO
     166!         ENDDO
     167!       else
    163168          teta(:,:)=phc(:,:)
    164        end if
     169!       end if
    165170     
    166171      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.10)) then
    … …  
    402407!====================================================================
    403408
     409      ztimestep=dt/real(ndt)
     410      do i=1,ndt
    404411
    405412!  Calcul de l,  km, au pas precedent
    406413      do k=2,nlay
    407                                                           do ig=1,ngrid
     414       do ig=1,ngrid
    408415!        print*,'SMML=',sm(ig,k),l(ig,k)
    409416         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
    … …  
    418425     s   (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
    419426! abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
    420          aa(ig,k)=aa1*dt/(delta(ig,k)*l(ig,k))
     427         aa(ig,k)=aa1*ztimestep/(delta(ig,k)*l(ig,k))
    421428!     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
    422429         qpre=sqrt(q2(ig,k))
    … …  
    440447
    441448!     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
    442                                                           enddo
     449       enddo
    443450      enddo
    444451
    … …  
    446453      q2(:,nlay+1)=q2(:,nlay)
    447454
     455      if (iflag_pbl .eq. 9) then
     456      do k=2,nlay
     457      do ig=1,ngrid
     458        zq=sqrt(q2(ig,k))
     459        km(ig,k)=l(ig,k)*zq*sm(ig,k)
     460        kn(ig,k)=km(ig,k)*alpha(ig,k)
     461        kq(ig,k)=l(ig,k)*zq*0.2
     462      enddo
     463      enddo
     464      km(:,nlay+1)=km(:,nlay)
     465      kn(:,nlay+1)=kn(:,nlay)
     466      kq(:,nlay+1)=kq(:,nlay)
     467
     468      q2(:,1)=q2(:,2)
     469       call vdif_q2(ztimestep,g,rconst,ngrid,nlay,plev,temp,kq,q2)
     470
     471      endif ! of if iflag_pbl eq 9
     472
     473      enddo !of i=1,ndt
     474
    448475      endif ! Fin du cas 8
    449476
    … …  
    453480!   Calcul des coefficients de melange
    454481!====================================================================
     482      if (iflag_pbl .ne. 9) then
    455483      do k=2,nlay
    456484!     print*,'k=',k
    … …  
    471499
    472500! Transport diffusif vertical de la TKE.
    473       if (iflag_pbl.ge.9) then
    474 !       print*,'YAMADA VDIF'
    475         q2(:,1)=q2(:,2)
    476         call vdif_q2(dt,g,rconst,ngrid,nlay,plev,temp,kq,q2)
     501!      if (iflag_pbl.ge.9) then
     502!!       print*,'YAMADA VDIF'
     503!        q2(:,1)=q2(:,2)
     504!        call vdif_q2(dt,g,rconst,ngrid,nlay,plev,temp,kq,q2)
     505!      endif
     506
    477507      endif
    478508
    … …  
    484514!   D'apres Holtslag Boville.
    485515
    486 ! MARS : deactivating pblhmin following F.H. concerns
    487 
    488 !                                                          do ig=1,ngrid
     516! MARS
     517!       callkmin=.true.
     518!       call getin("callkmin",callkmin)
     519!       IF (callkmin) THEN
     520                                                          do ig=1,ngrid
    489521!      coriol(ig)=1.e-4
    490522!      pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
    491 !                                                          enddo
    492 !
     523
     524       if (ngrid .eq. 1) then
     525       kminfact=0.3
     526       else
     527       kminfact=0.45
     528       endif
     529
     530       pblhmin(ig)=kminfact*0.07*MAX(ustar(ig),1.e-3)/1.e-4
     531                                                   enddo
    493532!      print*,'pblhmin ',pblhmin
    494533!CTest a remettre 21 11 02
    495534! test abd 13 05 02      if(0.eq.1) then
    496535!      if(0.eq.1) then
    497 !      do k=2,nlay
    498 !         do ig=1,ngrid
    499 !            if (teta(ig,2).gt.teta(ig,1)) then
    500 !               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
     536      do k=2,nlay
     537         do ig=1,ngrid
     538            if (teta(ig,2).gt.teta(ig,1)) then
     539               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
    501540!               kmin=kap*zlev(ig,k)*qmin
    502 !            else
    503 !               kmin=-1. ! kmin n'est utilise que pour les SL stables.
    504 !            endif
    505 !            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
     541               kmin=fl(zlev(ig,k),l0(ig),qmin**2,n2(ig,k))*qmin
     542            else
     543               kmin=-1. ! kmin n'est utilise que pour les SL stables.
     544            endif
     545            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
    506546!               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
    507547!     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
    … …  
    509549!               km(ig,k)=kmin
    510550!               kq(ig,k)=kmin
     551
     552               kn(ig,k)=kmin*alpha(ig,k)
     553               km(ig,k)=kmin
     554               kq(ig,k)=kmin*0.2
    511555!   la longueur de melange est suposee etre l= kap z
    512556!   K=l q Sm d'ou q2=(K/l Sm)**2
    513557!               q2(ig,k)=(qmin/sm(ig,k))**2
    514 !            endif
    515 !         enddo
    516 !      enddo
     558               q2(ig,k)=(kmin/
     559     &     (fl(zlev(ig,k),l0(ig),qmin**2,n2(ig,k))*sm(ig,k)))**2
     560            endif
     561         enddo
     562      enddo
    517563!      endif
     564
     565!      ENDIF
    518566
    519567!   Diagnostique pour stokage
Note: See TracChangeset for help on using the changeset viewer.