Changeset 2628 for trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90

Timestamp:

Feb 28, 2022, 6:46:07 PM (3 years ago)

Author:

abierjon

Message:

Mars GCM:
Big changes on mountain top dust flows for GCM6:

• the scheme now activates only in grid meshes that contain a mountain among a hard-written list, instead of every meshes. This is done to prevent strong artificial reinjections of dust in places that don't present huge converging slopes enabling the concentration of dust (ex: Valles Marineris, Hellas). Topdust is now also detrained as soon as it leaves the column it originated from.
• the list of the 19 allowed mountains is used by the subroutine topmons_setup in module topmons_mod, to compute a logical array contains_mons(ngrid). alpha_hmons and hsummit are also set up once and for all by this subroutine, which is called in physiq_mod's firstcall.
• contains_mons, alpha_hmons and hsummit are now saved variables of the module surfdat_h, and are called as such and not as arguments in the subroutines using them.
• the logical flag "slpwind" and the comments in the code have also been updated to the new terminology "mountain top dust flows", accordingly to ticket #71. The new flag read in callphys.def is "topflows".

AB

File:

• trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90 (modified) (17 diffs)

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

@@ -2 +2 @@
 
       IMPLICIT NONE
-
-!     sub-grid scale mountain mesh fraction
-      REAL, SAVE, ALLOCATABLE :: alpha_hmons(:)
-
-!$OMP THREADPRIVATE(alpha_hmons)
 
       CONTAINS
@@ -32 +27 @@
                                  clearsky,totcloudfrac,                &
 !             input sub-grid scale mountain
-                                 nohmons,hsummit,                      &
+                                 nohmons,                              &
 !             output
                                  pdqtop,wfin,dsodust,dsords,dsotop,    &
@@ -43 +38 @@
       USE dimradmars_mod, only: albedo,naerkind
       USE comsaison_h, only: dist_sol,mu0,fract
-      USE surfdat_h, only: emis,co2ice,hmons,summit
+      USE surfdat_h, only: emis,co2ice,hmons,summit,alpha_hmons, &
+                           hsummit,contains_mons
       USE callradite_mod, only: callradite
 
@@ -87 +83 @@
 !     input sub-grid scale mountain
       LOGICAL, INTENT(IN) :: nohmons
-      REAL, INTENT(IN) :: hsummit(ngrid)
 
 !--------------------------------------------------------
@@ -196 +191 @@
 
 !     Detrainment      
-      REAL coefdetrain(ngrid,nlayer)          ! coefficient for detrainment : % of stormdust detrained
+      REAL coefdetrain(ngrid,nlayer)          ! coefficient for detrainment : % of topdust detrained
       REAL dqdet_topdust_mass(ngrid,nlayer)   ! tendancy pdq topdust mass after detrainment only
       REAL dqdet_topdust_number(ngrid,nlayer) ! tendancy pdq topdust number after detrainment only
@@ -206 +201 @@
       ! **********************************************************************
       ! **********************************************************************
-      ! Parametrization of the entrainment by slope wind above the sub-grid
-      ! scale topography
+      ! Parametrization of the entrainment of dust by slope winds above the 
+      ! converging sub-grid scale topography ("mountain top dust flows")
       ! **********************************************************************
       ! **********************************************************************      
@@ -279 +274 @@
        ! 1.1. Call the second radiative transfer for topdust, obtain the extra heating
        ! *********************************************************************
+       
+       ! NB: since the only grid meshes that matter for the radiative transfer
+       !     are the ones that contain a mount (contains_mons(ig)=.true.),
+       !     it could be relevant to optimize the code by calling the RT
+       !     only for those grid meshes.
+       !     See Ticket #92 on trac.lmd.jussieu.fr/Planeto   
         CALL callradite(icount,ngrid,nlayer,nq,zday,zls,zq,albedo,     &
                  emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,   &
@@ -285 +286 @@
                  tau,aerosol,dsodust,tauscaling,dust_rad_adjust,    &
                  taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,rstormdust,rtopdust, &
-                 totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,&
+                 totstormfract,clearatm,dsords,dsotop,nohmons,&
                  clearsky,totcloudfrac)
        ! **********************************************************************
@@ -291 +292 @@
        ! **********************************************************************
         DO ig=1,ngrid
-          IF ( (mu0(ig) .gt. mu0lim) .and. (alpha_hmons(ig) .gt. 0.) ) THEN
+          IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN
+            !! mu0(ig)>mu0lim ensures that it is daytime
+            !! contains_mons=True ensures that there is a mount in the mesh and alpha_hmons>0
+          
             !! **********************************************************************
             !! Temperature profile above the mountain and in the close environment
@@ -353 +357 @@
               endif
             ENDDO
-          ENDIF ! IF ((mu0(ig) .gt. mu0lim) .and. alpha_hmons(ig) .gt. 0.)
+          ENDIF ! IF ((mu0(ig) .gt. mu0lim) .and. contains_mons(ig))
         ENDDO ! DO ig=1,ngrid
 
@@ -364 +368 @@
        ! **********************************************************************
        DO ig=1,ngrid
-         IF ( (mu0(ig) .gt. mu0lim) .and. (alpha_hmons(ig) .gt. 0.) ) THEN
+         IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN
            !! Positive buoyancy: negative vertical velocity entrains UP
            IF (dt_top(ig) .gt. 0.) THEN
@@ -477 +481 @@
 
         DO ig=1,ngrid
-          IF ( (mu0(ig) .gt. mu0lim) .and. (alpha_hmons(ig) .gt. 0.) ) THEN
+          IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN
             !! Total air mass within the PBL before entrainment (=> by PBL we mean between the surface and the layer where the vertical wind is maximum)
             masse_pbl(ig)=0.
@@ -489 +493 @@
        ! **********************************************************************
         DO ig=1,ngrid
-          IF ( (mu0(ig) .gt. mu0lim) .and. (alpha_hmons(ig) .gt. 0.) .and. (masse_pbl(ig) .gt. 0.) ) THEN
+          IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) .and. (masse_pbl(ig) .gt. 0.) ) THEN
             !! Transport of background dust + concentrated topdust above lwmax
             DO l=lwmax(ig),nlayer
@@ -583 +587 @@
           DO l=1,nlayer!-1
             !! Detrainment during the day
-            IF ( (mu0(ig) .gt. mu0lim) .and. (zq_topdust_mass(ig,l) .gt. zq_dust_mass(ig,l)*0.01)) THEN
+            IF ( (mu0(ig) .gt. mu0lim) .and. (zq_topdust_mass(ig,l) .gt. zq_dust_mass(ig,l)*0.01) .and. (contains_mons(ig)) ) THEN
                coefdetrain(ig,l)=1.*( rhobarz(ig,l+1)*abs(wfin(ig,l+1)) - rhobarz(ig,l)*abs(wfin(ig,l)) ) / masse(ig,l)
                !! Detrainment when abs(w(l)) > abs(w(l+1)), i.e. coefdetrain < 0
@@ -601 +605 @@
                !  dqdet_topdust_number(ig,l)=-(1.-exp(coefdetrain(ig,l)*ptimestep))*zq_dust_number(ig,l)/ptimestep
                endif
-            !! Full detrainment during the night imposed
+            !! Full detrainment imposed during the night or when topdust leaves its origin column (contains_mons=False)
             ELSE
                dqdet_topdust_mass(ig,l)=-zq_topdust_mass(ig,l)/ptimestep
@@ -681 +685 @@
                CALL WRITEDIAGFI(ngrid,'wfin_top', &
                 'wfin_top','',3,wfin(:,:))
+               CALL WRITEDIAGFI(ngrid,'hmons', &
+                'hmons','',2,hmons)
+               CALL WRITEDIAGFI(ngrid,'hsummit', &
+                'hsummit','',2,hsummit)
                CALL WRITEDIAGFI(ngrid,'alpha_hmons', &
                 'alpha_hmons','',2,alpha_hmons)
@@ -961 +969 @@
 
 
-      end subroutine van_leer
-
+      end subroutine van_leer   
+      
 !********************************************************************************
-       ! initialization module variables
-       subroutine ini_topmons_mod(ngrid,nlayer)
-       
-       implicit none
-       
-       integer, intent(in) :: ngrid
-       integer, intent(in) :: nlayer
-       
-       allocate(alpha_hmons(ngrid))
-
-       end subroutine ini_topmons_mod
-       
-       subroutine end_topmons_mod
-       
-       implicit none
-       
-       if (allocated(alpha_hmons)) deallocate(alpha_hmons)
-
-       end subroutine end_topmons_mod       
+      subroutine topmons_setup(ngrid)
+      ! Purpose:
+      ! 1) Fill the logical array contains_mons(:),
+      !    with contains_mons(ig)=True if there is
+      !    a mountain in the mesh ig
+      ! 2) Compute alpha_hmons(:) and hsummit(:)
+      use surfdat_h,only: phisfi,hmons,summit,base,&
+                          alpha_hmons,hsummit,contains_mons                                     
+      use comcstfi_h,only: pi,g
+      use planetwide_mod,only: planetwide_maxval, planetwide_minval
+      use geometry_mod,only: longitude_deg,latitude_deg,&
+                             boundslon,boundslat !boundslon/lat(ngrid,4) :
+                                                 !  |------------------------------|
+                                                 !  |north_west=2      north_east=1|
+                                                 !  |                              |
+                                                 !  |             (ig)             |
+                                                 !  |                              |
+                                                 !  |south_west=3      south_east=4| 
+                                                 !  |------------------------------| 
+                              
+      implicit none
+      integer,intent(in) :: ngrid
+      
+      ! Local variables
+      integer,parameter :: ntop_max = 19 ! total number of mounts written in the hmons list
+      real lon_top(ntop_max),lat_top(ntop_max) ! coordinates of the mounts (in deg)
+      ! Mountains list :
+      ! Olympus Mons,Ascraeus Mons,Elysium Mons,Arsia Mons,Pavonis Mons,
+      ! Hecates Tholus,Tharsis Tholus,Ceraunius Tholus,Alba Mons,Apollinaris Mons,
+      ! Albor Tholus,Biblis Tholus,Anseris Mons,Ulysses Tholus,Aeolis Mons,
+      ! Euripus Mons,Hadriacus Mons,Tyrrhenus Mons,Uranius Mons
+      !
+      ! NB: in 64x48 horiz. resolution, Biblis Tholus & Ulysses Tholus fall in the
+      !     same mesh, hence only Biblis Tholus is kept by the alpha_hmons computation
+      data lon_top/-134.,-104.5,146.9,-121.1,-113.4,&
+                    150.2,-90.8,-97.4,-109.6,174.4,&
+                    150.4,-124.6,86.6,-121.6,137.8,&
+                        105.,91.8,106.5,-92.2/
+      data lat_top/18.4,11.8,24.8,-8.4,-0.8,&
+                   31.8,13.4,24.,40.4,-9.3,&
+                   18.8,2.6,-29.8,2.9,-5.4,&
+                   -44.8,-32.1,-21.1,26.8/
+      integer,parameter :: ntop = 19 ! the topmons scheme is limited to the first ntop mounts
+      real :: boundslon_deg(4),boundslat_deg(4)
+      real :: hmax,hmin
+      integer :: ig,itop
+      
+      
+      
+      IF (ngrid.gt.1) THEN     
+        ! Sanity check
+        if (ntop.gt.ntop_max) then
+          call abort_physic("topmons_setup","Number of mountains ntop greater than ntop_max",1)
+        endif
+
+        ! Determine contains_mons
+        contains_mons(:)=.false.
+
+        do ig=1,ngrid
+          boundslon_deg(:)=boundslon(ig,:)/pi*180.
+          boundslat_deg(:)=boundslat(ig,:)/pi*180.
+
+          do itop=1,ntop
+            if ( (lon_top(itop).gt.boundslon_deg(2)).and.(lon_top(itop).lt.boundslon_deg(1)) &
+            .and.(lat_top(itop).gt.boundslat_deg(3)).and.(lat_top(itop).lt.boundslat_deg(2)) ) then
+              contains_mons(ig)=.true.
+              write(*,*) "topmons_setup: Found a mount at:"
+              write(*,*) "(",boundslon_deg(2),",",boundslat_deg(2),")  (",boundslon_deg(1),",",boundslat_deg(1),")"
+              write(*,*) "             ((",lon_top(itop),",",lat_top(itop),"))"
+              write(*,*) "(",boundslon_deg(3),",",boundslat_deg(3),")  (",boundslon_deg(4),",",boundslat_deg(4),")"
+            endif
+          enddo
+        enddo
+
+        ! Compute alpha_hmons
+        call planetwide_maxval(hmons,hmax)
+        call planetwide_minval(hmons,hmin)
+        do ig=1,ngrid
+          if (contains_mons(ig)) then
+          ! the mesh ig contains a mountain 
+            alpha_hmons(ig)= 0.5*(hmons(ig)-hmin)/(hmax-hmin)
+            ! Sanity check
+            if (alpha_hmons(ig).le.0) then
+              call abort_physic("topmons_setup","ERROR: alpha_hmons cannot be <0 "//&
+                                "if the mesh contains a mountain. Please check your "//&
+                                "formula or your start files.",1)
+            endif
+          else
+          ! the mesh ig doesn't contain a mountain
+            alpha_hmons(ig)= 0 
+          endif
+        enddo
+
+        ! Compute hsummit
+        hsummit(:)=summit(:)-phisfi(:)/g
+      
+      ELSE ! 1D case
+            hmin=0.
+            hmax=23162.1 !set here the height of the sub-grid scale topography
+            do ig=1,ngrid               
+              alpha_hmons(ig)= (hmons(ig)-hmin)/(hmax-hmin) !0.1*(hmons(ig)-hmin)/(hmax-hmin)
+              print*,"1D, hmons=",hmons(ig),"alpha=",alpha_hmons(ig)
+            enddo
+            
+            hsummit(:)=14000.
+      ENDIF ! (ngrid.gt.1)
+      
+      end subroutine topmons_setup
+!********************************************************************************   
       
       END MODULE topmons_mod