Changeset 544

Timestamp:

Feb 27, 2012, 10:44:32 AM (13 years ago)

Author:

acolaitis

Message:

27/02/12 == AC

Continuation of thermals setting, comparisons with mesoscale results on Case C
Added possibility to call gcm (or 1d) with constant prescribed sensible heat flux, in the spirit of direct comparisons with LES

... This is directly comparable to the variable tke_heat_flux in namelist.input
... Requires the use of yamada4.F from terrestrial GCM (mixes more, seem more numerically stable)
... Usually requires high timesteps (>1000) to avoids crashes. Best approach is to compare
height of first model level z1 and teta_1 between LES and 1D, and increase the timesteps until results
between the two models are comparable (might require a slitghly different tke_heat_flux between the two models
due to difference in vertical diffusion schemes and subgrid effects)
... Syntax for use is to add "tke_heat_flux = ..." in callphys.def

Corrected some stuff with tke transport in thermals (which is anyway deactivated but one day maybe...)

Location:

trunk/LMDZ.MARS

Files:

• README (modified) (1 diff)
• libf/phymars/callkeys.h (modified) (2 diffs)
• libf/phymars/calltherm_interface.F90 (modified) (3 diffs)
• libf/phymars/inifis.F (modified) (1 diff)
• libf/phymars/physiq.F (modified) (4 diffs)
• libf/phymars/thermcell_main_mars.F90 (modified) (8 diffs)
• libf/phymars/vdifc.F (modified) (4 diffs)
• libf/phymars/yamada4.F (added)

trunk/LMDZ.MARS/README

@@ -1414 +1414 @@
 >> Code about 10% faster [tests carried out with mesoscale summer Tharsis run]
 >> Check with a GCM run that results similar
+
+== 27/02/12 == AC
+>> Continuation of thermals setting, comparisons with mesoscale results on Case C
+>> Added possibility to call gcm (or 1d) with constant prescribed sensible heat flux, in the spirit of direct comparisons with LES
+   ... This is directly comparable to the variable tke_heat_flux in namelist.input
+   ... Requires the use of yamada4.F from terrestrial GCM (mixes more, seem more numerically stable)
+   ... Usually requires high timesteps (>1000) to avoids crashes. Best approach is to compare
+   height of first model level z1 and teta_1 between LES and 1D, and increase the timesteps until results 
+   between the two models are comparable (might require a slitghly different tke_heat_flux between the two models
+   due to difference in vertical diffusion schemes and subgrid effects)
+   ... Syntax for use is to add "tke_heat_flux = ..." in callphys.def
+>> Corrected some stuff with tke transport in thermals (which is anyway deactivated but one day maybe...)

trunk/LMDZ.MARS/libf/phymars/callkeys.h

@@ -18 +18 @@
      &   ,dustbin,nqchem_min,nltemodel,nircorr
      
-      COMMON/callkeys_r/topdustref,solarcondate,semi,alphan
+      COMMON/callkeys_r/topdustref,solarcondate,semi,alphan,            &
+     &   tke_heat_flux
      
       LOGICAL callrad,calldifv,calladj,callcond,callsoil,               &
@@ -36 +37 @@
       real alphan
       real solarcondate
+      real tke_heat_flux
 
       integer iddist

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

@@ -158 +158 @@
        endif
 
-!       dtke_thermals=.false. !! default setting
-!       !call getin("dtke_thermals",dtke_thermals)
-!
-!       IF(dtke_thermals) THEN
-!          DO l=1,nlayermx
-!              q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1))
-!          ENDDO
-!       ENDIF
+       dtke_thermals=.false. !! default setting
+       call getin("dtke_thermals",dtke_thermals)
+       IF(dtke_thermals) THEN
+          DO l=1,nlayermx
+              q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1))
+          ENDDO
+       ENDIF
 
 ! **********************************************************************
@@ -338 +337 @@
       endif
 
+      if (dtke_thermals) then
+      modname='tke'
+      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
+     &     ,fm_therm,entr_therm,detr_therm,  &
+     &    masse,q2_therm,dq2_therm,modname,zdz)
+      endif
+
       DO ig=1,ngridmx
          hfmax(ig)=MAXVAL(heatFlux(ig,:)+heatFlux_down(ig,:))
@@ -369 +375 @@
            endif
 
-!           IF(dtke_thermals) THEN
-!              DO l=2,nlayermx
-!                 pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))
-!              ENDDO
-!  
-!              pbl_dtke(:,1)=0.5*dq2_therm(:,1)
-!              pbl_dtke(:,nlayermx+1)=0.
-!           ENDIF
+           IF(dtke_thermals) THEN
+              DO l=2,nlayermx
+                 pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))
+              ENDDO
+  
+              pbl_dtke(:,1)=0.5*dq2_therm(:,1)
+              pbl_dtke(:,nlayermx+1)=0.
+           ENDIF
 
 

trunk/LMDZ.MARS/libf/phymars/inifis.F

@@ -203 +203 @@
          write(*,*) " topdustref = ",topdustref
 
+         write(*,*) "Prescribed surface thermal flux (H/(rho*cp),K m/s)"
+         tke_heat_flux=0. ! default value
+         call getin("tke_heat_flux",tke_heat_flux)
+         write(*,*) " tke_heat_flux = ",tke_heat_flux
+         write(*,*) " 0 means the usual schemes are computing"
+
          write(*,*) "call radiative transfer ?"
          callrad=.true. ! default value

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

@@ -11 +11 @@
 #endif
      $            )
+
 
       IMPLICIT NONE
@@ -338 +339 @@
 
 c Variables for PBL
-
+      REAL zz1(ngridmx)
       REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
       REAL, SAVE :: wstar(ngridmx)
@@ -769 +770 @@
 #endif
 
+         IF (tke_heat_flux .ne. 0.) THEN
+             zz1(:)=(pt(:,1)+pdt(:,1)*ptimestep)*(r/g)*
+     &                 (-alog(pplay(:,1)/pplev(:,1)))
+             pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
+         ENDIF
 
 c        Calling vdif (Martian version WITH CO2 condensation)
@@ -1997 +2003 @@
          call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
          call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
+!         call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",              &
+!     &              "K.s-1",1,dtrad/zpopsk)
+!        call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
+!     &                   'w.m-2',1,zdtsw/zpopsk)
+!        call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
+!     &                   'w.m-2',1,zdtlw/zpopsk)
+
 ! or output in diagfi.nc (for testphys1d)
          call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)

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

@@ -113 +113 @@
       REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
       LOGICAL activecell(ngridmx),activetmp(ngridmx)
-      REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega
+      REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega,adalim
       INTEGER tic
 
@@ -381 +381 @@
       b1inv=b1
       omega=0.
+      adalim=0.
 
 ! One good config for 34/35 levels
@@ -389 +390 @@
 ! Best configuration for 222 levels:
 
-      omega=0.06
-      b1=0.
-      a1=1.
-      a1inv=0.25*a1
-      b1inv=0.0002
+!      omega=0.06
+!      b1=0.
+!      a1=1.
+!      a1inv=0.25*a1
+!      b1inv=0.0002
+!!
+!!
+!!      ae=0.9*ae
 
 ! Best config for norad 222 levels:
-
-!       omega=0.06
-!       a1=1.
-!       b1=0.
-!       a1inv=a1
-!       be=1.1*be
-!       ad = 0.0004
-!       b1inv=0.00035
+! with yamada4 and alim at sqrt(zlev)
+
+       omega=0.06
+       a1=1.
+       b1=0.
+       a1inv=a1
+       be=1.1*be
+       ad = 0.0004
+       b1inv=0.00035
+       adalim=0.
+
+      b1inv=0.00025
 
 ! Trying stuff :
+
+!      omega=0.04
+!!      b1=0.
+!      a1=1.
+!      a1inv=a1
+!      b1inv=0.0005689
+!!      be=1.1*be
+!!      ae=0.96*ae
+
 
 !       omega=0.06
@@ -431 +448 @@
             if (ztv(ig,1)>=(ztv(ig,2))) then
                alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.)  &
-     &                       *sqrt(zlev(ig,2))
+!     &                       *sqrt(zlev(ig,2))
 !     &                       /sqrt(zlev(ig,2))
-!      &                       *zlev(ig,2)
+      &                       *zlev(ig,2)
                lalim(ig)=2
                alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,1)
@@ -444 +461 @@
            if (ztv(ig,l)>(ztv(ig,l+1)+0.) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.)) then ! .and. (zlev(ig,l+1) .lt. 1000.)) then
                alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
-     &                       *sqrt(zlev(ig,l+1))
-!      &                       *zlev(ig,2)
+!     &                       *sqrt(zlev(ig,l+1))
+!     &                       /sqrt(zlev(ig,l+1))
+      &                       *zlev(ig,l+1)
                 lalim(ig)=l+1
                alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
@@ -562 +580 @@
           if(zbuoy(ig,l) .gt. 0.) then
              if(l .lt. lalim(ig)) then
-                detr_star(ig,l)=0.
+
+!                detr_star(ig,l)=0.
+                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
+            &  adalim
              else
 
@@ -628 +649 @@
 !---------------------------------------------------------------------------
 
-      DO tic=0,0  ! internal convergence loop
+      DO tic=0,5  ! internal convergence loop
       activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10
       do ig=1,ngridmx
@@ -674 +695 @@
           if(zbuoy(ig,l) .gt. 0.) then
              if(l .lt. lalim(ig)) then
-                detr_star(ig,l)=0.
+
+!                detr_star(ig,l)=0.
+                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
+            &  adalim
+
              else
                  detr_star(ig,l) = f_star(ig,l)*zdz*              &

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

@@ -10 +10 @@
 #endif
      &                )
-       USE ioipsl_getincom
       IMPLICIT NONE
 
@@ -369 +368 @@
 c    ** schema de diffusion turbulente dans la couche limite
 c       ---------------------------------------------------- 
+       IF (tke_heat_flux .eq. 0.) THEN
 
        CALL vdif_kc(ptimestep,g,pzlev,pzlay
@@ -374 +374 @@
      &              ,pq2,zkv,zkh,zq)
 
-!      pt(:,:)=ph(:,:)*ppopsk(:,:)
-!      CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
-!     s   ,pzlev,pzlay,pu,pv,ph,zcdv_true,pq2,zkv,zkh,zkq,ust
-!     s   ,8)
+      ELSE
+
+      pt(:,:)=zh(:,:)*ppopsk(:,:)
+      CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
+     s   ,pzlev,pzlay,pu,pv,ph,zcdv_true,pq2,zkv,zkh,zkq,ust
+     s   ,8)
+
+      ENDIF
 
       if ((doubleq).and.(ngrid.eq.1)) then
@@ -523 +527 @@
 
       z4st=4.*5.67e-8*ptimestep
+      IF (tke_heat_flux .eq. 0.) THEN
       DO ig=1,ngrid
          zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
       ENDDO
+      ELSE
+         zdplanck(:)=0.
+      ENDIF
 
 ! calcul de zc et zd pour la couche top en prenant en compte le terme