Changeset 621 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Apr 16, 2012, 12:00:40 PM (13 years ago)
Author:
acolaitis
Message:

Thermals. Finally corrected advection scheme... The bug was coming from a non conservation of fluxes, not from the advection scheme itself.

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r619 r621  
    6565      REAL lmax_real(ngridmx)
    6666      REAL masse(ngridmx,nlayermx)
    67       REAL zdz(ngridmx,nlayermx)
    6867      LOGICAL qtransport_thermals,dtke_thermals
    6968      INTEGER l,ig,iq,ii(1)
    70       CHARACTER (LEN=20) :: modname
    7169
    7270!--------------------------------------------------------
    … …  
    112110! **********************************************************************
    113111
    114       lmax(:)=0.
     112      lmax(:)=0
    115113      pdu_th(:,:)=0.
    116114      pdv_th(:,:)=0.
    … …  
    236234         zzw2(:,:)=0.
    237235         zmax(:)=0.
    238          lmax(:)=0.
     236         lmax(:)=0
    239237         d_t_the(:,:)=0.
    240238         d_u_the(:,:)=0.
    … …  
    311309!****************************************************************
    312310
     311      lmax(:)=nint(lmax_real(:))
     312
    313313! Now that we have computed total entrainment and detrainment, we can
    314314! advect u, v, and q in thermals. (theta already advected). We can do
    … …  
    321321      enddo
    322322
    323 ! thickness of layers
    324       do l=1,nlayermx
    325          zdz(:,l)=zzlev(:,l+1)-zzlev(:,l)
    326       enddo
    327 
    328       modname='momentum'
    329323      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
    330      &      ,fm_therm,entr_therm,zfraca,  &
    331      &     masse,zu,d_u_ajs,modname,zdz)
     324     &      ,fm_therm,entr_therm,  &
     325     &     masse,zu,d_u_ajs)
    332326
    333327      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
    334      &       ,fm_therm,entr_therm,zfraca,  &
    335      &     masse,zv,d_v_ajs,modname,zdz)
     328     &       ,fm_therm,entr_therm,  &
     329     &     masse,zv,d_v_ajs)
    336330
    337331      if (nqmx .ne. 0.) then
    338       modname='tracer'
    339332      DO iq=1,nqmx
    340333      if (iq .ne. ico2) then
    341334      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
    342      &     ,fm_therm,entr_therm,zfraca,  &
    343      &    masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),modname,zdz)
     335     &     ,fm_therm,entr_therm,  &
     336     &    masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq))
    344337      endif
    345338      ENDDO
    … …  
    347340
    348341      if (dtke_thermals) then
    349       modname='tke'
    350342      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
    351      &     ,fm_therm,entr_therm,zfraca,  &
    352      &    masse,q2_therm,dq2_therm,modname,zdz)
     343     &     ,fm_therm,entr_therm,  &
     344     &    masse,q2_therm,dq2_therm)
    353345      endif
    354346
    … …  
    357349         wmax(ig)=MAXVAL(zw2(ig,:))
    358350      ENDDO
    359 
    360       lmax(:)=nint(lmax_real(:))
    361351
    362352! **********************************************************************

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

    r620 r621  
    11      subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm0,entr0,  &
    2      &    alpha,masse0,q_therm,dq_therm,charvar,zdz)
     2     &    masse0,q_therm,dq_therm)
    33      implicit none
    44
    … …  
    2323      REAL, INTENT(IN) ::entr0(ngridmx,nlayermx)
    2424      REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx)
    25       CHARACTER (LEN=20), INTENT(IN) :: charvar
    2625      REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)
    27       REAL, INTENT(IN) :: zdz(ngridmx,nlayermx)
    28       REAL, INTENT(IN) :: alpha(ngridmx,nlayermx+1)
    2926
    3027! ============================ OUTPUTS ===========================
    … …  
    3633      REAL q(ngridmx,nlayermx)
    3734      REAL qa(ngridmx,nlayermx)
    38       REAL qd(ngridmx,nlayermx)
    3935      INTEGER ig,k
    40       REAL gammaf(ngridmx,nlayermx),fmid(ngridmx,nlayermx)
    41       REAL gammae(ngridmx,nlayermx)
    42       REAL alphamid(ngridmx,nlayermx)
    43       REAL Z(2,2,ngridmx,nlayermx)
    44       REAL X(2,1,ngridmx,nlayermx)
     36      REAL detr(ngridmx,nlayermx),entr(ngridmx,nlayermx)
     37      REAL wqd(ngridmx,nlayermx+1)
     38      REAL zzm
     39
    4540! =========== Init ==============================================
    4641
    4742      qa(:,:)=q_therm(:,:)
    4843      q(:,:)=q_therm(:,:)
     44      entr(:,:)=entr0(:,:)
     45      detr(:,:)=0.
    4946
    50       do ig=1,ngridmx
    51          do k=1, nlayermx
    52             if (fm0(ig,k)+entr0(ig,k) .gt. 0.) then
    53               gammaf(ig,k)=fm0(ig,k)/(fm0(ig,k)+entr0(ig,k))
    54             else
    55               gammaf(ig,k)=0.
    56             endif
    57             if (fm0(ig,k+1) .gt. 0.) then
    58               gammae(ig,k)=(fm0(ig,k)+entr0(ig,k))/fm0(ig,k+1)
    59             else
    60               gammae(ig,k)=1.
     47      do k=1,nlayermx
     48         do ig=1,ngridmx
     49            detr(ig,k)=fm0(ig,k)-fm0(ig,k+1)+entr(ig,k)
     50            if (detr(ig,k).lt.0.) then
     51               entr(ig,k)=entr(ig,k)-detr(ig,k)
     52               detr(ig,k)=0.
    6153            endif
    6254         enddo
    6355      enddo
    64       do k=1,nlayermx
    65         do ig=1, ngridmx
    66           Z(1,1,ig,k)=1.
    67           Z(1,2,ig,k)=0.
    68           Z(2,1,ig,k)=0.
    69           Z(2,2,ig,k)=1.         
    70         enddo
    71       enddo
    72       do k=2,nlayermx
    73         do ig=1, ngridmx
    74           if (gammaf(ig,k) .gt. 0.) then
    75           Z(1,1,ig,k)=1./gammae(ig,k-1) + (gammaf(ig,k)-1.)*            &
    76      &               (1.-gammae(ig,k-1))/(gammae(ig,k-1)*gammaf(ig,k))
    77           Z(1,2,ig,k)=-(1.-gammae(ig,k-1))                              &
    78      &                        /(gammae(ig,k-1)*gammaf(ig,k))
    79           Z(2,1,ig,k)=(gammaf(ig,k)-1.)/gammaf(ig,k)
    80           Z(2,2,ig,k)=1./gammaf(ig,k)
    81           else
    82           Z(1,1,ig,k)=1./gammae(ig,k-1) -1.
    83           Z(1,2,ig,k)=1.
    84           Z(2,:,ig,k)=0.
    85           endif
    8656
    87         enddo
     57! =========== Updraft ============================================
     58
     59      do ig=1,ngrid
     60         qa(ig,1)=q(ig,1)
    8861      enddo
    8962
    90       X(1,1,:,:)=q_therm(:,:)
    91       X(2,1,:,:)=q_therm(:,:)
    92 
    93       do k=nlayermx,2,-1
    94         do ig=1, ngridmx
    95        X(:,1,ig,k-1)=MATMUL(Z(:,:,ig,k),X(:,1,ig,k))
    96 
    97         enddo
     63      do k=2,nlayermx
     64         do ig=1,ngridmx
     65            if ((fm0(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
     66     &         1.e-5*masse0(ig,k)) then
     67         qa(ig,k)=(fm0(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
     68     &     /(fm0(ig,k+1)+detr(ig,k))
     69            else
     70               qa(ig,k)=q(ig,k)
     71            endif
     72         enddo
    9873      enddo
    9974
    100       q(:,:)=X(1,1,:,:)
    101       qa(:,:)=X(2,1,:,:)
     75! =========== Subsidence ==========================================
     76
     77      do k=2,nlayermx
     78         do ig=1,ngridmx
     79
     80! Schema avec advection sur plus qu'une maille.
     81            zzm=masse0(ig,k)/ptimestep
     82            if (fm0(ig,k)>zzm) then
     83               wqd(ig,k)=(zzm*q(ig,k)+(fm0(ig,k)-zzm)*q(ig,k+1))
     84            else
     85               wqd(ig,k)=fm0(ig,k)*q(ig,k)
     86            endif
     87         enddo
     88      enddo
     89      do ig=1,ngrid
     90         wqd(ig,1)=0.
     91         wqd(ig,nlayermx+1)=0.
     92      enddo
     93
     94
     95! ====== dq ======================================================
    10296
    10397      dq_therm(:,:)=0.
    104       fmid(:,:)=0.
    105       alphamid(:,:)=0.
    10698
    107         do ig=1, ngridmx
    108            do k=1,nlayermx
    109               fmid(ig,k) = 0.5*(fm0(ig,k)+fm0(ig,k+1))
    110               alphamid(ig,k) = 0.5*(alpha(ig,k)+alpha(ig,k+1))
    111            enddo
    112         enddo
     99         do k=1,nlayermx-1
     100           q(:,k)=q(:,k)+         &
     101     & (detr(:,k)*qa(:,k)-entr(:,k)*q(:,k) &
     102     &       -wqd(:,k)+wqd(:,k+1))  &
     103     &               *ptimestep/masse0(:,k)
     104         enddo
    113105
    114         do ig=1, ngridmx
    115            do k=1,nlayermx-1
    116               dq_therm(ig,k)=-(1./masse0(ig,k))*(  &
    117      &   (1.-alphamid(ig,k+1))*fmid(ig,k+1)*(qa(ig,k+1)-q(ig,k+1)) -   &
    118      &   (1.-alphamid(ig,k))*fmid(ig,k)*(qa(ig,k)-q(ig,k))          ) &
    119      &       /zdz(ig,k)
     106         do k=1,nlayermx-1
     107          dq_therm(:,k)=(q(:,k)-q_therm(:,k))/ptimestep
     108         enddo
    120109
    121            enddo
    122         enddo
    123110      return
    124111      end

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

    r619 r621  
    175175      data firstcall/.true./
    176176      REAL zhc(ngridmx,nlayermx)
    177       REAL zdzfull(ngridmx,nlayermx)
    178177      REAL ratiom(ngridmx,nlayermx)
    179178
    … …  
    411410       ad = 0.0004
    412411       adalim=0.
     412
    413413       b1inv=0.00025
     414!       b1inv=b1
     415
    414416!       b1inv = 0.0003
    415417
    … …  
    12281230      zdthladj(:,:)=0.
    12291231
     1232      if (1 .eq. 0) then
     1233      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
     1234     &     ,fm,entr,  &
     1235     &    masse,ztv,zdthladj)
     1236      else
     1237
     1238
    12301239      do ig=1,ngridmx
    12311240         if(lmax(ig) .gt. 1) then
    … …  
    12441253         endif
    12451254      enddo
     1255
     1256      endif
    12461257
    12471258! ------------------------------------------------------------------
    … …  
    15281539      if (ico2.ne.0) then
    15291540!      if (nqmx .ne. 0.) then
    1530       do l=1,nlayermx
    1531          zdzfull(:,l)=zlev(:,l+1)-zlev(:,l)
    1532       enddo
    1533 
    1534       modname='tracer'
    15351541      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
    1536      &     ,fm,entr,fraca,  &
    1537      &    masse,pq(:,:,ico2),pdqadj(:,:,ico2),modname,zdzfull)
     1542     &     ,fm,entr,  &
     1543     &    masse,pq(:,:,ico2),pdqadj(:,:,ico2))
    15381544!      endif
    15391545
Note: See TracChangeset for help on using the changeset viewer.