Changeset 2594 for LMDZ5/branches/testing/libf/phylmd/grid_noro_m.F90

Timestamp:

Jul 18, 2016, 9:41:10 PM (8 years ago)

Author:

Laurent Fairhead

Message:

Merged trunk changes r2545:2589 into testing branch

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phylmd/grid_noro_m.F90 (modified) (3 diffs)

LMDZ5/branches/testing/libf/phylmd/grid_noro_m.F90

@@ -3 +3 @@
 !*******************************************************************************
 
+  USE print_control_mod, ONLY: lunout
+  USE assert_eq_m,       ONLY: assert_eq
   PRIVATE
-  PUBLIC :: grid_noro
+  PUBLIC :: grid_noro, grid_noro0
 
 
@@ -45 +47 @@
 !    on the edge are contained in input cells.
 !===============================================================================
-  USE assert_eq_m, ONLY: assert_eq
-  USE print_control_mod, ONLY: lunout
-  IMPLICIT NONE 
-  REAL, PARAMETER :: epsfra = 1.e-5
+  IMPLICIT NONE
 !-------------------------------------------------------------------------------
 ! Arguments:
@@ -315 +314 @@
 !-------------------------------------------------------------------------------
 !
+SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
+!
+!===============================================================================
+! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
+!          without any call to physics subroutines.
+!===============================================================================
+  IMPLICIT NONE 
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
+  REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
+  REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
+  REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
+  REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="grid_noro0"
+  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
+  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
+  REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
+  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
+  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
+  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
+  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
+  LOGICAL :: masque_lu
+  INTEGER :: i, ii, imdp, imar, iext
+  INTEGER :: j, jj, jmdp, jmar, nn
+  REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
+  REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
+!-------------------------------------------------------------------------------
+  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
+  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
+  imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
+  jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
+! IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
+! IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
+  iext=imdp/10
+  xpi = ACOS(-1.)
+  rad = 6371229.
+
+!--- ARE WE USING A READ MASK ?
+  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
+  WRITE(lunout,*)'Masque lu: ',masque_lu
+
+!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
+  ALLOCATE(xusn(imdp+2*iext))
+  xusn(1     +iext:imdp  +iext)=xd(:)
+  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
+  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
+
+  ALLOCATE(yusn(jmdp+2))
+  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
+  yusn(2:jmdp+1)=yd(:)
+  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
+
+  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
+  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
+  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
+  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
+  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
+  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
+  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
+  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
+
+!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
+  ALLOCATE(a(imar+1),b(imar+1))
+  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
+  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
+  a(1)=x(1)-(x(2)-x(1))/2.0
+  a(2:imar+1)= b(1:imar)
+
+  ALLOCATE(c(jmar),d(jmar))
+  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
+  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
+  c(1)=y(1)-(y(2)-y(1))/2.0
+  c(2:jmar)=d(1:jmar-1)
+
+!--- INITIALIZATIONS:
+  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
+  ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
+
+!--- SUMMATION OVER GRIDPOINT AREA
+  zleny=xpi/REAL(jmdp)*rad
+  xincr=xpi/REAL(jmdp)/2.
+  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
+  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
+  DO ii = 1, imar+1
+    DO jj = 1, jmar
+      DO j = 2,jmdp+1 
+        zlenx  =zleny  *COS(yusn(j))
+        zbordnor=(xincr+c(jj)-yusn(j))*rad
+        zbordsud=(xincr-d(jj)+yusn(j))*rad
+        weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
+        IF(weighy/=0) THEN
+          DO i = 2, imdp+2*iext-1
+            zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
+            zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
+            weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
+            IF(weighx/=0)THEN
+              num_tot(ii,jj)=num_tot(ii,jj)+1.0
+              IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
+              weight(ii,jj)=weight(ii,jj)+weighx*weighy
+              zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
+            END IF
+          END DO
+        END IF
+      END DO
+    END DO
+  END DO
+
+!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
+  IF(.NOT.masque_lu) THEN
+    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
+  END IF
+  nn=COUNT(weight(:,1:jmar-1)==0.0)
+  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
+  WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
+
+!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
+  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
+  WHERE(mask>=0.1) mask_tmp = 1.
+  WHERE(weight(:,:)/=0.0)
+    zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
+    zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
+  END WHERE
+  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
+
+!--- Values at poles
+  zphi(imar+1,:)=zphi(1,:)
+
+  zweinor=SUM(weight(1:imar,   1),DIM=1)
+  zweisud=SUM(weight(1:imar,jmar),DIM=1)
+  zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
+  zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
+  zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
+
+END SUBROUTINE grid_noro0
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
 SUBROUTINE MVA9(x)
 !