Changeset 284 for trunk/LMDZ.MARS/libf/phymars

Timestamp:

Sep 7, 2011, 3:20:37 PM (13 years ago)

Author:

acolaitis

Message:

--- AC 07/09/2011 ---

• Added new flag for the Richardson-based surface layer :

callrichsl, can be changed in callphys.def

One should always use the thermals model when using this surface layer model.
Somes cases (weakly unstable with low winds), when not using thermals, won't be well represented by the
Richardson surface layer. This stands for Mesoscale and Gcm but not for LES model.

Correct configs :

callrichsl = .true.
calltherm = .true.

callrichsl = .false.
calltherm = .false.

callrichsl = .false.
calltherm = .true.

Previously unstable config :

callrichsl = .true.
calltherm = .false.

• To be able to run without thermals and with the new surface layer, a modification has been made to

physiq.F to account for gustiness in GCM and MESOSCALE for negative Richardson, so that :

callrichsl = .true.
calltherm = .false.

can now be used without problems, but is not recommended.

• Consequently, callrichsl = .false. is now the default configuration for thermals.

We recall the available options in callphys.def for thermals :

outptherm = BOOLEAN (.false. by default) : outputs thermals related quantities (lots of diagfi)
nsplit_thermals = INTEGER (50 by default in gcm, 2 in mesoscale) : subtimestep for thermals model.

It is recommended to use at least 40 in the gcm, and at least 2 in the mesoscale.
The user can lower these values but should check it's log for anomalies or errors regarding
tracer transport in the thermals, or "granulosity" in the outputs for wmax, lmax and hfmax.

---

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• callkeys.h (modified) (2 diffs)
• calltherm_interface.F90 (modified) (3 diffs)
• calltherm_mars.F90 (modified) (1 diff)
• inifis.F (modified) (1 diff)
• meso_inc/meso_inc_les.F (modified) (1 diff)
• physiq.F (modified) (3 diffs)
• vdif_cd.F (modified) (4 diffs)
• vdifc.F (modified) (4 diffs)

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

@@ -12 +12 @@
      &   ,callg2d,linear,rayleigh,tracer,active,doubleq,submicron       &
      &   ,lifting,callddevil,scavenging,sedimentation,activice,water    &
-     &   ,caps,photochem,calltherm,outptherm,callslope
+     &   ,caps,photochem,calltherm,outptherm,callrichsl,callslope
      
       COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
@@ -24 +24 @@
      &   ,callnirco2,callnlte,callthermos,callconduct,                  &
      &    calleuv,callmolvis,callmoldiff,thermochem,thermoswater        &
-     &   ,calltherm,outptherm,callslope
+     &   ,calltherm,outptherm,callrichsl,callslope
 
 

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

@@ -71 +71 @@
       REAL buoyancyEst(ngridmx,nlayermx)
       REAL hfmax(ngridmx),wmax(ngridmx)
-
-      REAL tstart,tstop
 
 !---------------------------------------------------------
@@ -135 +133 @@
          endif
 
-         call cpu_time(tstart)
-
          CALL calltherm_mars(ptimestep,zzlev,zzlay &
      &      ,pplay,pplev,pphi &
@@ -146 +142 @@
      &      ,zpopsk,ztla,heatFlux,heatFlux_down &
      &      ,buoyancyOut,buoyancyEst,hfmax,wmax)
-          call cpu_time(tstop)
-         print*,'TOTAL elapsed time in thermals : ',tstop-tstart
-
 
 ! Accumulation des  tendances. On n'accumule pas les quantités de traceurs car celle ci n'a pas du changer

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

@@ -99 +99 @@
          r_aspect_thermals=0.7
 #else
-         nsplit_thermals=40
+         nsplit_thermals=50
          r_aspect_thermals=2.
 #endif

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

@@ -252 +252 @@
           stop
          endif
+
+         write(*,*) "call Richardson-based surface layer ?"
+         callrichsl=.false. ! default value
+         call getin("callrichsl",callrichsl)
+         write(*,*) " callrichsl = ",callrichsl
 
          write(*,*) "call CO2 condensation ?"

trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F

@@ -1 @@
-         DO ig=1,ngrid
-          !! sensible heat flux in W/m2
-!          hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig)
+         if (callrichsl .eq. .false.) then
+       
+             DO ig=1,ngrid
+!! sensible heat flux in W/m2
+
+             hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig)
+
+!! u star in similarity theory in m/s
+             ust(ig) = 0.4
+     .               * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) )
+     .               / log( 1.E+0 + zzlay(ig,1)/z0_default )
+             ENDDO
+
+         else
+
+            DO ig=1,ngrid
 
 ! New SL parametrization, correct formulation for hfx :
 
-          hfx(ig) = (pplay(ig,1)/(r*pt(ig,1)))*cpp
-     &    *sqrt((pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)))
-     &    *zcdh(ig)*(tsurf(ig)-zh(ig,1))
+            hfx(ig) = (pplay(ig,1)/(r*pt(ig,1)))*cpp
+     &        *sqrt(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
+     &        + (1.0*wmax_th(ig))**2)
+     &        *zcdh(ig)*(tsurf(ig)-zh(ig,1))
 
-          !! u star in similarity theory in m/s
-!          ust(ig) = 0.4
-!     .               * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) )
-!     .               / log( 1.E+0 + zzlay(ig,1)/z0_default )
 
 ! New SL parametrization, ust is more accurately computed in vdif_cd :
-        ust(ig) = sqrt(zcdv(ig)*(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)))
+            ust(ig) = sqrt(zcdv(ig)*
+     &   (pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + (1.0*wmax_th(ig))**2)
+     &                     )
 
-         ENDDO   
+            ENDDO   
+
+         endif  !of if callrichsl
+
 !         write (*,*) 'PHYS HFX cp zdts', hfx(100), zflubid(100), 
 !     .       capcal(100), 

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

@@ -694 +694 @@
          enddo
 
+c Treatment of a special case : using new surface layer (Richardson based)
+c without using the thermals in gcm and mesoscale can yield problems in
+c weakly unstable situations when winds are near to 0. For those cases, we add
+c a unit subgrid gustiness. Remember that thermals should be used we using the 
+c Richardson based surface layer model.
+
+#ifdef MESOSCALE
+      IF (flag_LES .eq. .false.) THEN
+        IF ((calltherm .eq. .false.) .and. (callrichsl .eq. .true.)) THEN
+          DO ig=1, ngridmx
+             IF (zh(ig,1) .lt. tsurf(ig)) THEN
+               wmax_th(ig)=1.
+             ENDIF        
+          ENDDO
+        ENDIF
+      ENDIF
+#else
+      IF ((calltherm .eq. .false.) .and. (callrichsl .eq. .true.)) THEN
+        DO ig=1, ngridmx
+          IF (zh(ig,1) .lt. tsurf(ig)) THEN
+            wmax_th(ig)=1.
+          ENDIF
+        ENDDO
+      ENDIF    
+#endif
+
 c        Calling vdif (Martian version WITH CO2 condensation)
          CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
@@ -1433 +1459 @@
         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
      &                  co2ice)
+
 c         call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
 c     &                  "K",2,zt(1,7))
-         call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
-     &                  fluxsurf_lw)
-         call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
-     &                  fluxsurf_sw_tot)
-         call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
-     &                  fluxtop_lw)
-         call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
-     &                  fluxtop_sw_tot)
+c         call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
+c     &                  fluxsurf_lw)
+c         call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
+c     &                  fluxsurf_sw_tot)
+c         call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
+c     &                  fluxtop_lw)
+c         call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
+c     &                  fluxtop_sw_tot)
          call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
         call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
@@ -1734 +1761 @@
          endif
 
-! ---
-
-
-
          if(calltherm) then
 

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

@@ -36 +36 @@
 
 #include "comcstfi.h"
+#include "callkeys.h"
 
 c   Arguments:
@@ -87 +88 @@
 c   ------------------
 
-      reynolds(:)=0.
-
 c Original formulation :
 
-c      DO ig=1,ngrid
-c         z1=1.E+0 + pz(ig,1)/pz0(ig)
-c         zcd0=karman/log(z1)
-c         zcd0=zcd0*zcd0
-c         pcdv(ig)=zcd0
-c         pcdh(ig)=zcd0
-c      ENDDO
+      if(callrichsl .eq. .false.) then
+
+      DO ig=1,ngrid
+         z1=1.E+0 + pz(ig,1)/pz0(ig)
+         zcd0=karman/log(z1)
+         zcd0=zcd0*zcd0
+         pcdv(ig)=zcd0
+         pcdh(ig)=zcd0
+      ENDDO
      
 !      print*,'old : cd,ch; ',pcdv,pcdh
+      else
+
+      reynolds(:)=0.
 
 c New formulation (AC) :
@@ -224 +228 @@
 ! Some useful diagnostics :
 
-!        call WRITEDIAGFI(ngrid,'RiB',
+!       call WRITEDIAGFI(ngrid,'RiB',
 !     &              'Bulk Richardson nb','no units',
 !     &                         2,rib)
@@ -241 +245 @@
 
 
+      endif !of if call richardson surface layer
+
       RETURN
       END

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

@@ -86 +86 @@
       REAL zc(ngridmx,nlayermx),zd(ngridmx,nlayermx)
       REAL zcst1
-      REAL zu2
+      REAL zu2(ngridmx)
 
       EXTERNAL SSUM,SCOPY
@@ -234 +234 @@
      &             ,zcdv_true,zcdh_true)
 
-      DO ig=1,ngrid
-
-        zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
-     &             +(1.*wmax(ig))**2        ! subgrid gustiness added by quadratic interpolation with a coeff beta, here beta=1. (tuned from
+
+        IF (callrichsl .eq. .true.) then
+           zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)
+     &             +(1.*wmax(:))**2        ! subgrid gustiness added by quadratic interpolation with a coeff beta, here beta=1. (tuned from
                                             ! LES comparisons. This parameter is linked to the thermals settings)
-        
-        zcdv(ig)=zcdv_true(ig)*sqrt(zu2)     ! 1 / bulk aerodynamic momentum conductance 
-        zcdh(ig)=zcdh_true(ig)*sqrt(zu2)     ! 1 / bulk aerodynamic heat conductance
+        ELSE
+           zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)
+        ENDIF        
+
+        zcdv(:)=zcdv_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance 
+        zcdh(:)=zcdh_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
                                              ! these are the quantities to be looked at when comparing surface layers of different models
-      ENDDO
 
 ! Some usefull diagnostics for the new surface layer parametrization :
@@ -254 +256 @@
 !     &                         2,zcdh_true)
 !        call WRITEDIAGFI(ngridmx,'u*',
-!    &              'friction velocity','m/s',
-!     &                         2,sqrt(zcdv_true(:))*sqrt(zu2))
+!     &              'friction velocity','m/s',
+!     &                         2,sqrt(zcdv_true(:))*sqrt(zu2(:)))
 !        call WRITEDIAGFI(ngridmx,'T*',
 !     &              'friction temperature','K',
@@ -479 +481 @@
 c     ----------------------------------------------------------------
         DO ig=1,ngrid
-          zu2=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1)
-          zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
-          zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
+          zu2(ig)=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1)
+          zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig))
+          zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig))
         ENDDO