Changeset 3468

Timestamp:

Mar 14, 2019, 5:43:22 PM (5 years ago)

Author:

Laurent Fairhead

Message:

Modifications so that point close to the equator (and the height of the stratosphere)
can be found on an unstructured grid
YM

File:

• LMDZ6/trunk/libf/phylmd/orografi_strato.F90 (modified) (3 diffs)

LMDZ6/trunk/libf/phylmd/orografi_strato.F90

@@ -1819 +1819 @@
   USE mod_phys_lmdz_para
   USE mod_grid_phy_lmdz
+  USE geometry_mod
   IMPLICIT NONE
 
@@ -1832 +1833 @@
   INTEGER jk
   REAL zpr, ztop, zsigt, zpm1r
-  REAL :: pplay_glo(klon_glo, nlev)
-  REAL :: paprs_glo(klon_glo, nlev+1)
+  INTEGER :: cell,ij,nstra_tmp,nktopg_tmp
+  REAL :: current_dist, dist_min,dist_min_glo
 
   ! *       1.    SET THE VALUES OF THE PARAMETERS
@@ -1848 +1849 @@
   ! old  ZSIGT=0.85
 
-  CALL gather(pplay, pplay_glo)
-  CALL bcast(pplay_glo)
-  CALL gather(paprs, paprs_glo)
-  CALL bcast(paprs_glo)
-
-  DO jk = 1, nlev
-    zpm1r = pplay_glo(klon_glo/2+1, jk)/paprs_glo(klon_glo/2+1, 1)
-    IF (zpm1r>=zsigt) THEN
-      nktopg = jk
-    END IF
-    zpm1r = pplay_glo(klon_glo/2+1, jk)/paprs_glo(klon_glo/2+1, 1)
-    IF (zpm1r>=ztop) THEN
-      nstra = jk
-    END IF
-  END DO
-
+
+!ym Take the point at equator close to (0,0) coordinates.
+  dist_min=360
+  cell=-1
+  DO ij=1,klon
+    current_dist=sqrt(longitude_deg(ij)**2+latitude_deg(ij)**2)
+    current_dist=current_dist*(1+(1e-10*ind_cell_glo(ij))/klon_glo) ! For point unicity
+    IF (dist_min>current_dist) THEN
+      dist_min=current_dist
+      cell=ij    
+    ENDIF  
+  ENDDO
+  
+  CALL reduce_min(dist_min,dist_min_glo)
+  IF (dist_min/=dist_min_glo) cell=-1
+!ym in future find the point at equator close to (0,0) coordinates.
+
+  nktopg_tmp=nktopg
+  nstra_tmp=nstra
+  
+  IF (cell/=-1) THEN
+
+    DO jk = 1, nlev
+      zpm1r = pplay(cell+1, jk)/paprs(cell+1, 1)
+      IF (zpm1r>=zsigt) THEN
+        nktopg_tmp = jk
+      END IF
+      IF (zpm1r>=ztop) THEN
+        nstra_tmp = jk
+      END IF
+    END DO
+  ELSE
+    nktopg_tmp=0
+    nstra_tmp=0
+  ENDIF
+  
+  CALL reduce_sum(nktopg_tmp,nktopg)
+  CALL bcast(nktopg)
+  CALL reduce_sum(nstra_tmp,nstra)
+  CALL bcast(nstra)
+  
   ! inversion car dans orodrag on compte les niveaux a l'envers
   nktopg = nlev - nktopg + 1