Changeset 2930 for trunk/LMDZ.MARS/libf

Timestamp:

Apr 2, 2023, 9:08:46 PM (20 months ago)

Author:

emillour

Message:

Mars PCM:
Adapt topmons_setup to handle both cartesian and icosahedral grids
EM

File:

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

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

@@ -989 +989 @@
       use comcstfi_h,only: pi,g
       use planetwide_mod,only: planetwide_maxval, planetwide_minval
+      use mod_grid_phy_lmdz, only: nvertex !number of encompassing vertex
       use geometry_mod,only: longitude_deg,latitude_deg,&
-                             boundslon,boundslat !boundslon/lat(ngrid,4) :
+                             boundslon,boundslat ! if nvertex==4 then we have
+                                                 ! boundslon/lat(ngrid,4) :
                                                  !  |------------------------------|
                                                  !  |north_west=2      north_east=1|
@@ -998 +1000 @@
                                                  !  |south_west=3      south_east=4| 
                                                  !  |------------------------------| 
+                                                 ! if nvertex==6 (isocahedral grid)
+                                                 ! there is no dedicated ordering
                               
       implicit none
@@ -1022 +1026 @@
                    -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 :: boundslon_deg(nvertex),boundslat_deg(nvertex)
       real :: hmax,hmin
+      real :: minlon,maxlon ! min/max encompassing vertex longitude (deg)
+      real :: minlat,maxlat ! min/max encompassing vertex latitude (deg)
       integer :: ig,itop
       
@@ -1040 +1046 @@
           boundslon_deg(:)=boundslon(ig,:)/pi*180.
           boundslat_deg(:)=boundslat(ig,:)/pi*180.
+          
+          minlon=minval(boundslon_deg)
+          maxlon=maxval(boundslon_deg)
+          minlat=minval(boundslat_deg)
+          maxlat=maxval(boundslat_deg)
 
           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
+            if ((lon_top(itop).ge.minlon).and.(lon_top(itop).le.maxlon).and. &
+                (lat_top(itop).ge.minlat).and.(lat_top(itop).le.maxlat)) 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),")"
+              write(*,*) "topmons_setup: Found a mount for ig:",ig
+              write(*,*) " at : lon=",longitude_deg(ig)," lat=",latitude_deg(ig)
+              if (nvertex==4) then
+               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
             endif
-          enddo
-        enddo
+          enddo ! of do itop=1,ntop
+        enddo ! of do ig=1,ngrid
 
         ! Compute alpha_hmons
         call planetwide_maxval(hmons,hmax)
         call planetwide_minval(hmons,hmin)
+
+        ! Add a sanity check in case there is no topography:
+        if (hmin==hmax) then
+          write(*,*) "topmons_setup: hmin=",hmin," hmax=",hmax
+          write(*,*) "topmons_setup: thus no topflows !!!"
+          contains_mons(:)=.false.
+        endif
+
         do ig=1,ngrid
           if (contains_mons(ig)) then
@@ -1070 +1092 @@
             alpha_hmons(ig)= 0 
           endif
-        enddo
+        enddo ! of do ig=1,ngrid
 
         ! Compute hsummit