source: LMDZ5/trunk/libf/dynphy_lonlat/grid_atob_m.f90 @ 4268

Last change on this file since 4268 was 2532, checked in by dcugnet, 8 years ago

Bug fix: For high resolutions or particular zoom center locations, an output cell (LMDZ grid) with no matching input cells ("landiceref.nc" file grid) caused a crash.
File size: 10.8 KB
RevLine 
[2293]1!*******************************************************************************
2!
3MODULE grid_atob_m
4!
5!*******************************************************************************
6
7  USE assert_eq_m, ONLY: assert_eq
8
9  PRIVATE
10  PUBLIC :: grille_m, rugosite, sea_ice, rugsoro
11
12CONTAINS
13
14!-------------------------------------------------------------------------------
15!
16SUBROUTINE fine2coarse(x_i, y_i, x_o, y_o, d_o1, d_i, msk, d_o2)
17!
18!-------------------------------------------------------------------------------
19  IMPLICIT NONE
20!-------------------------------------------------------------------------------
21! Arguments:
[2299]22  REAL,              INTENT(IN)  :: x_i(:), y_i(:) !-- INPUT  X&Y COOR. (mi)(ni)
[2361]23  REAL,              INTENT(IN)  :: x_o(:), y_o(:) !-- OUTPUT X&Y COOR. (mo)(no)
[2299]24  REAL,              INTENT(OUT) :: d_o1(:,:)      !-- OUTPUT FIELD     (mo,no)
25  REAL,    OPTIONAL, INTENT(IN)  :: d_i (:,:)      !-- INPUT FIELD      (mi,ni)
[2361]26  LOGICAL, OPTIONAL, INTENT(IN)  :: msk (:,:)      !-- MASK             (mo,no)
[2299]27  REAL,    OPTIONAL, INTENT(OUT) :: d_o2(:,:)      !-- OUTPUT FOR d_i^2 (mo,no)
[2293]28!-------------------------------------------------------------------------------
29! Local variables:
30  CHARACTER(LEN=256) :: modname="fine2coarse"
[2532]31  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1,m2, n1,n2, mx,my, nn, i,j
32  LOGICAL :: li, lo, first=.TRUE.
33  REAL    :: inc, cpa, spa, crlo(SIZE(x_i))
34  REAL, SAVE :: pi, hpi
[2361]35  INTEGER, DIMENSION(SIZE(x_o),SIZE(y_o)) :: num_tot
36  LOGICAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: found, mask
[2532]37  REAL,    DIMENSION(SIZE(x_i),SIZE(y_i)) :: dist
[2361]38  REAL,    DIMENSION(SIZE(x_o))           :: a, b
39  REAL,    DIMENSION(SIZE(y_o))           :: c, d
40  REAL,    PARAMETER :: thresh=1.E-5
[2293]41!-------------------------------------------------------------------------------
[2532]42  IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF
[2361]43  mi=SIZE(x_i); ni=SIZE(y_i); mo=SIZE(x_o); no=SIZE(y_o)
44  m1=m1; m2=mo; mx=mo; IF(PRESENT(msk)) mx=SIZE(msk,1)
45  n1=ni; n2=no; my=no; IF(PRESENT(msk)) my=SIZE(msk,2)
[2293]46  li=PRESENT(d_i ); IF(li) THEN; m1=SIZE(d_i ,1); n1=SIZE(d_i ,2); END IF
47  lo=PRESENT(d_o2); IF(lo) THEN; m2=SIZE(d_o2,1); n2=SIZE(d_o2,2); END IF
[2361]48  mi=assert_eq(mi,m1,TRIM(modname)//" mi")
49  ni=assert_eq(ni,n1,TRIM(modname)//" ni")
50  mo=assert_eq(mo,m2,mx,SIZE(d_o1,1),TRIM(modname)//" mo")
51  no=assert_eq(no,n2,my,SIZE(d_o1,2),TRIM(modname)//" no")
52  mask(:,:)=.TRUE.; IF(PRESENT(msk)) mask(:,:)=msk(:,:)
[2293]53
54!--- COMPUTE CELLS INTERFACES COORDINATES OF OUTPUT GRID
55  b(mo)=x_o(mo)+(x_o(mo)-x_o(mo-1))/2.0; b(1:mo-1)=(x_o(1:mo-1)+x_o(2:mo))/2.0
56  d(no)=y_o(no)+(y_o(no)-y_o(no-1))/2.0; d(1:no-1)=(y_o(1:no-1)+y_o(2:no))/2.0
57  a(1 )=x_o(1 )-(x_o(2 )-x_o(1   ))/2.0; a(2:mo  )=   b(1:mo-1)
58  c(1 )=y_o(1 )-(y_o(2 )-y_o(1   ))/2.0; c(2:no  )=   d(1:no-1)
59
60!--- ACCUMULATE INPUT POINTS ON OUTPUT GRID
[2299]61  d_o1(:,:)=0.; num_tot(:,:)=0; inc=1.0
[2293]62  IF(lo) d_o2(:,:)=0.
63  DO ji = 1, ni
64    DO ii = 1, mi
[2299]65      IF(li) inc=d_i(ii,ji)
[2293]66      DO jo = 1, no
67        IF((y_i(ji)-c(jo)<thresh.OR.y_i(ji)-d(jo)>thresh).AND.   &
68           (y_i(ji)-c(jo)>thresh.OR.y_i(ji)-d(jo)<thresh)) CYCLE
69        DO io = 1, mo
70          IF((x_i(ii)-a(io)<thresh.OR.x_i(ii)-b(io)>thresh).AND. &
71             (x_i(ii)-a(io)>thresh.OR.x_i(ii)-b(io)<thresh)) CYCLE
72          num_tot(io,jo)=num_tot(io,jo)+1
[2361]73          IF(mask(io,jo)) d_o1(io,jo)=d_o1(io,jo)+inc
[2293]74          IF(.NOT.lo) CYCLE
[2361]75          IF(mask(io,jo)) d_o2(io,jo)=d_o2(io,jo)+inc*inc
[2293]76        END DO
77      END DO
78    END DO
79  END DO
80
81!--- CHECK INPUT POINTS HAVE BEEN FOUND IN EACH OUTPUT CELL
82  found(:,:)=num_tot(:,:)/=0
[2299]83  WHERE(found.AND.mask) d_o1(:,:)=d_o1(:,:)/REAL(num_tot(:,:))
[2293]84  IF(PRESENT(d_o2)) THEN
[2299]85    WHERE(found.AND.mask) d_o2(:,:)=d_o2(:,:)/REAL(num_tot(:,:))
[2293]86    RETURN
87  END IF
[2294]88  nn=COUNT(found); IF(nn==0) RETURN
[2293]89
90!--- MISSING POINTS ; USE DISTANCE ON THE SPHERE TO FIND NEAREST POINT nr(2)
91  DO io = 1, mo
92    DO jo = 1, no
93      IF(found(io,jo)) CYCLE
94!      IF(prt_level>=1) PRINT*, "Problem: point out of domain (i,j)=", io,jo
[2532]95      crlo(:)=COS(x_o(io)-x_i(:))     !--- COS of points 1 and 2 angle
96      cpa=COS(y_o(jo)); spa=SIN(y_o(jo))
97      DO j=1,ni; dist(:,j)=ACOS(spa*SIN(y_i(j))+cpa*COS(y_i(j))*crlo(:)); END DO
[2293]98      nr=MINLOC(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr
[2299]99      inc=1.0; IF(li) inc=d_i(nr(1),nr(2))
[2361]100      IF(mask(io,jo)) d_o1(io,jo)=inc
[2293]101    END DO
102  END DO
103
104END SUBROUTINE fine2coarse
105!
106!-------------------------------------------------------------------------------
107
108
109!-------------------------------------------------------------------------------
110!
111SUBROUTINE grille_m(xdata, ydata, entree, x, y, sortie)
112!
113!-------------------------------------------------------------------------------
114! Author:  Z.X. Li (april 1st 1994) (see also A. Harzallah and L. Fairhead)
115!-------------------------------------------------------------------------------
116! Purpose: Naive method to regrid on a coarser grid.
117!   Value at a new point is an average of the old points lcoated in a zone
118!   surrounding the considered new point.
119!   No ponderation is used (see grille_p)
120!
121!           (c)
122!        ----d-----
123!        | . . . .|
124!        |        |
125!     (b)a . * . .b(a)
126!        |        |
127!        | . . . .|
128!        ----c-----
129!           (d)
130!
131!-------------------------------------------------------------------------------
132  IMPLICIT NONE
133!-------------------------------------------------------------------------------
134! Arguments:
135  REAL, INTENT(IN)  :: xdata(:), ydata(:)       !--- INPUT  FIELD X AND Y COORD.
136  REAL, INTENT(IN)  :: entree(SIZE(xdata),SIZE(ydata)) !--- INPUT  FIELD
137  REAL, INTENT(IN)  :: x(:), y(:)               !--- OUTPUT FIELD X AND Y COORD.
138  REAL, INTENT(OUT) :: sortie(SIZE(x),SIZE(y))  !--- OUTPUT FIELD
139!-------------------------------------------------------------------------------
[2299]140  CALL fine2coarse(xdata,ydata,x,y,sortie,entree)
[2293]141
142END SUBROUTINE grille_m
143!
144!-------------------------------------------------------------------------------
145
146
147!-------------------------------------------------------------------------------
148!
149SUBROUTINE rugosite(xdata, ydata, entree, x, y, sortie, mask)
150!
151!-------------------------------------------------------------------------------
152! Author:  Z.X. Li (april 1st 1994)
153!-------------------------------------------------------------------------------
[2299]154! Purpose: Remap rugosity length ; constant value (0.001) on oceans.
[2293]155! Naive method  (see grille_m)
156!-------------------------------------------------------------------------------
157  IMPLICIT NONE
158!-------------------------------------------------------------------------------
159! Arguments:
160  REAL, INTENT(IN)  :: xdata(:), ydata(:)      !--- INPUT  FIELD X AND Y COORD.
161  REAL, INTENT(IN)  :: entree(SIZE(xdata),SIZE(ydata)) !--- INPUT  FIELD
162  REAL, INTENT(IN)  :: x(:), y(:)              !--- OUTPUT FIELD X AND Y COORD.
163  REAL, INTENT(OUT) :: sortie(SIZE(x),SIZE(y)) !--- OUTPUT FIELD
164  REAL, INTENT(IN)  :: mask  (SIZE(x),SIZE(y)) !--- MASK
165!-------------------------------------------------------------------------------
[2299]166  CALL fine2coarse(xdata,ydata,x,y,sortie,LOG(entree))
[2293]167  WHERE(NINT(mask)==1)
[2299]168    sortie(:,:)=EXP(sortie(:,:))
[2293]169  ELSE WHERE
[2299]170    sortie(:,:)=0.001
[2293]171  END WHERE
172
173END SUBROUTINE rugosite
174!
175!-------------------------------------------------------------------------------
176
177
178!-------------------------------------------------------------------------------
179!
180SUBROUTINE sea_ice(xdata, ydata, glace01, x, y, frac_ice)
181!
182!-------------------------------------------------------------------------------
183! Author:  Z.X. Li (april 1st 1994)
184! Purpose: Regrid ice indicator (0 or 1) on a coarser grid to get an ice fract.
185! field (between 0 and 1).
186! Naive method  (see grille_m)
187!-------------------------------------------------------------------------------
188  IMPLICIT NONE
189!-------------------------------------------------------------------------------
190! Arguments:
191  REAL, INTENT(IN)  :: xdata(:), ydata(:)      !--- INPUT  FIELD X AND Y COORD.
192  REAL, INTENT(IN)  :: glace01(:,:)            !--- INPUT  FIELD
193  REAL, INTENT(IN)  :: x(:), y(:)              !--- OUTPUT FIELD X AND Y COORD.
194  REAL, INTENT(OUT) :: frac_ice(SIZE(x),SIZE(y)) !--- OUTPUT FIELD
195!-------------------------------------------------------------------------------
[2299]196  CALL fine2coarse(xdata,ydata,x,y,frac_ice,msk=NINT(glace01)==1)
[2293]197
198END SUBROUTINE sea_ice
199!
200!-------------------------------------------------------------------------------
201
202
203!-------------------------------------------------------------------------------
204!
205SUBROUTINE rugsoro(xrel, yrel, relief, xmod, ymod, rugs)
206!
207!-------------------------------------------------------------------------------
208! Author:  Z.X. Li (april 1st 1994) ; Modifications: august 23rd 1995.
209! Purpose: Compute rugosity due to orography by using standard dev in a 1x1 cell
210!-------------------------------------------------------------------------------
211  IMPLICIT NONE
212!-------------------------------------------------------------------------------
213! Arguments:
214  REAL, INTENT(IN)  :: xrel(:), yrel(:)        !--- INPUT  FIELD X AND Y COORD.
215  REAL, INTENT(IN)  :: relief(:,:)             !--- INPUT  FIELD
216  REAL, INTENT(IN)  :: xmod(:), ymod(:)        !--- OUTPUT FIELD X AND Y COORD.
217  REAL, INTENT(OUT) :: rugs(SIZE(xmod),SIZE(ymod)) !--- OUTPUT FIELD
218!-------------------------------------------------------------------------------
219! Local variable:
220  INTEGER           :: k, nn
221  INTEGER, PARAMETER:: itmp=360, jtmp=180
[2299]222  REAL  :: out(SIZE(xmod),SIZE(xmod)), amin, amax
223  REAL  :: cham1tmp(itmp,jtmp), cham2tmp(itmp,jtmp), xtmp(itmp), ytmp(jtmp)
[2293]224!-------------------------------------------------------------------------------
225
226!--- TEST INPUT FILE FITS FOR ONE DEGREE RESOLUTION
227  xtmp(:)=4.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(itmp),k=1,itmp)]
228  ytmp(:)=2.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(jtmp),k=1,jtmp)]
229  CALL fine2coarse(xrel,yrel,xtmp,ytmp,cham1tmp,relief,d_o2=cham2tmp)
230  cham2tmp(:,:)=cham2tmp(:,:)-cham1tmp(:,:)**2
231  nn=COUNT(cham2tmp<=-7.5)
232  IF(nn/=0) THEN
233    PRINT*,'Problem for rugsoro ; std**2 < -7.5 for several points: ',nn
234    CALL ABORT_GCM("", "", 1)
235  END IF
236  nn=COUNT(cham2tmp<0)
237  IF(nn/=0) PRINT*,'Problem for rugsoro ; std**2 < 0. for several points: ',nn
[2299]238  WHERE(cham2tmp<0.0) cham2tmp=0.0
[2293]239  cham2tmp(:,:)=SQRT(cham2tmp(:,:))
240  amin=MINVAL(cham2tmp); amax=MAXVAL(cham2tmp)
241  PRINT*, 'Ecart-type 1x1:', amin, amax
242
243!--- COMPUTE RUGOSITY AT REQUIRED SCALE
[2299]244  WHERE(cham2tmp<0.001) cham2tmp=0.001
[2293]245  CALL fine2coarse(xtmp,ytmp,xmod,ymod,out,REAL(LOG(cham2tmp)))
246  out=EXP(out)
247  amin=MINVAL(out); amax=MAXVAL(out)
248  PRINT*, 'Ecart-type du modele:', amin, amax
[2299]249  out=out/amax*20.0
[2293]250  amin=MINVAL(out); amax=MAXVAL(out)
251  PRINT*, 'Longueur de rugosite du modele:', amin, amax
252  rugs=REAL(out)
253
254END SUBROUTINE rugsoro
255!
256!-------------------------------------------------------------------------------
257
258END MODULE grid_atob_m
259!
260!*******************************************************************************
261
Note: See TracBrowser for help on using the repository browser.