source: LMDZ5/trunk/libf/dynlonlat_phylonlat/grid_atob_m.F90 @ 2294

Last change on this file since 2294 was 2294, checked in by fhourdin, 9 years ago

Introduction d'un seuil sur la probabilité de non déclenchement
stochastique random_notrig_max
Si le tirage aleatoire uniforme entre 0 et 1 est > random_notrig_max
on ne déclenche pas.
random_notrig_max=1-epsilon avec epsilon petit.

Corrections dans grid_atob_m.F90 pour la compilation Linux.
A vérifier.
File size: 12.6 KB
Line 
1!*******************************************************************************
2!
3MODULE grid_atob_m
4!
5!*******************************************************************************
6
7  USE assert_eq_m, ONLY: assert_eq
8  REAL, SAVE :: pi, deg2rad
9
10  PRIVATE
11  PUBLIC :: grille_m, rugosite, sea_ice, rugsoro
12
13CONTAINS
14
15!-------------------------------------------------------------------------------
16!
17SUBROUTINE fine2coarse(x_i, y_i, x_o, y_o, d_o1, d_i, msk, d_o2)
18!
19!-------------------------------------------------------------------------------
20  IMPLICIT NONE
21!-------------------------------------------------------------------------------
22! Arguments:
23  REAL,                       INTENT(IN)  :: x_i(:), y_i(:)  !- IN  X&Y COORD.
24  REAL,                       INTENT(IN)  :: x_o(:), y_o(:)  !- OUT X&Y COORD.
25  DOUBLE PRECISION,           INTENT(OUT) :: d_o1(:,:)       !- OUT FLD (mo,no)
26  REAL,             OPTIONAL, INTENT(IN)  :: d_i (:,:)       !- INP FLD (mi,ni)
27  LOGICAL,          OPTIONAL, INTENT(IN)  :: msk (:,:)       !- MASK    (mi,ni)
28  DOUBLE PRECISION, OPTIONAL, INTENT(OUT) :: d_o2(:,:)       !- OUT FOR d_i^2
29!-------------------------------------------------------------------------------
30! Local variables:
31  CHARACTER(LEN=256) :: modname="fine2coarse"
32  DOUBLE PRECISION   :: inc
33  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1, n1, m2, n2, nn
34  INTEGER :: num_tot(SIZE(x_o),SIZE(y_o))
35  LOGICAL :: found(SIZE(x_o),SIZE(y_o)), li
36  LOGICAL :: mask (SIZE(x_i),SIZE(y_i)), lo
37  REAL    :: dist (SIZE(x_o),SIZE(y_o))
38  REAL    :: a(SIZE(x_o)), b(SIZE(x_o)), c(SIZE(y_o)), d(SIZE(y_o))
39  REAL, PARAMETER :: thresh=1.E-5
40!-------------------------------------------------------------------------------
41  mask(:,:)=.TRUE.; IF(PRESENT(msk)) mask(:,:)=msk(:,:)
42  mi=SIZE(x_i); m1=mi; ni=SIZE(y_i); n1=ni
43  mo=SIZE(x_o); m2=mo; no=SIZE(y_o); n2=no
44  li=PRESENT(d_i ); IF(li) THEN; m1=SIZE(d_i ,1); n1=SIZE(d_i ,2); END IF
45  lo=PRESENT(d_o2); IF(lo) THEN; m2=SIZE(d_o2,1); n2=SIZE(d_o2,2); END IF
46  mi=assert_eq(mi,m1,SIZE(mask,1),TRIM(modname)//" mi")
47  ni=assert_eq(ni,n1,SIZE(mask,2),TRIM(modname)//" ni")
48  mo=assert_eq(mo,m2,SIZE(d_o1,1),TRIM(modname)//" mo")
49  no=assert_eq(no,n2,SIZE(d_o1,2),TRIM(modname)//" no")
50
51!--- COMPUTE CELLS INTERFACES COORDINATES OF OUTPUT GRID
52  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
53  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
54  a(1 )=x_o(1 )-(x_o(2 )-x_o(1   ))/2.0; a(2:mo  )=   b(1:mo-1)
55  c(1 )=y_o(1 )-(y_o(2 )-y_o(1   ))/2.0; c(2:no  )=   d(1:no-1)
56
57!--- ACCUMULATE INPUT POINTS ON OUTPUT GRID
58  d_o1(:,:)=0.; num_tot(:,:)=0; inc=1.0d0
59  IF(lo) d_o2(:,:)=0.
60  DO ji = 1, ni
61    DO ii = 1, mi
62      IF(li) inc=DBLE(d_i(ii,ji))
63      DO jo = 1, no
64        IF((y_i(ji)-c(jo)<thresh.OR.y_i(ji)-d(jo)>thresh).AND.   &
65           (y_i(ji)-c(jo)>thresh.OR.y_i(ji)-d(jo)<thresh)) CYCLE
66        DO io = 1, mo
67          IF((x_i(ii)-a(io)<thresh.OR.x_i(ii)-b(io)>thresh).AND. &
68             (x_i(ii)-a(io)>thresh.OR.x_i(ii)-b(io)<thresh)) CYCLE
69          num_tot(io,jo)=num_tot(io,jo)+1
70          IF(mask(ii,ji)) d_o1(io,jo)=d_o1(io,jo)+inc
71          IF(.NOT.lo) CYCLE
72          IF(mask(ii,ji)) d_o2(io,jo)=d_o2(io,jo)+inc*inc
73        END DO
74      END DO
75    END DO
76  END DO
77
78!--- CHECK INPUT POINTS HAVE BEEN FOUND IN EACH OUTPUT CELL
79  found(:,:)=num_tot(:,:)/=0
80  WHERE(found.AND.mask) d_o1(:,:)=d_o1(:,:)/DBLE(num_tot(:,:))
81  IF(PRESENT(d_o2)) THEN
82    WHERE(found.AND.mask) d_o2(:,:)=d_o2(:,:)/DBLE(num_tot(:,:))
83    RETURN
84  END IF
85  nn=COUNT(found); IF(nn==0) RETURN
86
87!--- MISSING POINTS ; USE DISTANCE ON THE SPHERE TO FIND NEAREST POINT nr(2)
88  DO io = 1, mo
89    DO jo = 1, no
90      IF(found(io,jo)) CYCLE
91!      IF(prt_level>=1) PRINT*, "Problem: point out of domain (i,j)=", io,jo
92      CALL dist_sphe(x_o(io),y_o(jo),x_i,y_i,dist(:,:))
93      nr=MINLOC(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr
94      inc=1.0; IF(li) inc=DBLE(d_i(nr(1),nr(2)))
95      IF(mask(nr(1),nr(2))) d_o1(io,jo)=inc
96    END DO
97  END DO
98
99END SUBROUTINE fine2coarse
100!
101!-------------------------------------------------------------------------------
102
103
104!-------------------------------------------------------------------------------
105!
106SUBROUTINE grille_m(xdata, ydata, entree, x, y, sortie)
107!
108!-------------------------------------------------------------------------------
109! Author:  Z.X. Li (april 1st 1994) (see also A. Harzallah and L. Fairhead)
110!-------------------------------------------------------------------------------
111! Purpose: Naive method to regrid on a coarser grid.
112!   Value at a new point is an average of the old points lcoated in a zone
113!   surrounding the considered new point.
114!   No ponderation is used (see grille_p)
115!
116!           (c)
117!        ----d-----
118!        | . . . .|
119!        |        |
120!     (b)a . * . .b(a)
121!        |        |
122!        | . . . .|
123!        ----c-----
124!           (d)
125!
126!-------------------------------------------------------------------------------
127  IMPLICIT NONE
128!-------------------------------------------------------------------------------
129! Arguments:
130  REAL, INTENT(IN)  :: xdata(:), ydata(:)       !--- INPUT  FIELD X AND Y COORD.
131  REAL, INTENT(IN)  :: entree(SIZE(xdata),SIZE(ydata)) !--- INPUT  FIELD
132  REAL, INTENT(IN)  :: x(:), y(:)               !--- OUTPUT FIELD X AND Y COORD.
133  REAL, INTENT(OUT) :: sortie(SIZE(x),SIZE(y))  !--- OUTPUT FIELD
134!-------------------------------------------------------------------------------
135! Local variable:
136  DOUBLE PRECISION  :: out(SIZE(x),SIZE(y))
137!-------------------------------------------------------------------------------
138!  CALL fine2coarse(xdata,ydata,x,y,out,DBLE(entree))
139  CALL fine2coarse(xdata,ydata,x,y,out,entree)
140  sortie=REAL(out)
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!-------------------------------------------------------------------------------
154! Purpose: From topography field, compute ocean/land mask (land: 1 ; ocean: 0)
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!-------------------------------------------------------------------------------
166! Local variable:
167  DOUBLE PRECISION  :: out   (SIZE(x),SIZE(y))
168!-------------------------------------------------------------------------------
169  CALL fine2coarse(xdata,ydata,x,y,out,LOG(entree))
170  WHERE(NINT(mask)==1)
171    out(:,:)=EXP(out(:,:))
172  ELSE WHERE
173    out(:,:)=0.001
174  END WHERE
175  sortie=REAL(out)
176
177END SUBROUTINE rugosite
178!
179!-------------------------------------------------------------------------------
180
181
182!-------------------------------------------------------------------------------
183!
184SUBROUTINE sea_ice(xdata, ydata, glace01, x, y, frac_ice)
185!
186!-------------------------------------------------------------------------------
187! Author:  Z.X. Li (april 1st 1994)
188! Purpose: Regrid ice indicator (0 or 1) on a coarser grid to get an ice fract.
189! field (between 0 and 1).
190! Naive method  (see grille_m)
191!-------------------------------------------------------------------------------
192  IMPLICIT NONE
193!-------------------------------------------------------------------------------
194! Arguments:
195  REAL, INTENT(IN)  :: xdata(:), ydata(:)      !--- INPUT  FIELD X AND Y COORD.
196  REAL, INTENT(IN)  :: glace01(:,:)            !--- INPUT  FIELD
197  REAL, INTENT(IN)  :: x(:), y(:)              !--- OUTPUT FIELD X AND Y COORD.
198  REAL, INTENT(OUT) :: frac_ice(SIZE(x),SIZE(y)) !--- OUTPUT FIELD
199!-------------------------------------------------------------------------------
200! Local variable:
201  DOUBLE PRECISION  :: out     (SIZE(x),SIZE(y))
202!-------------------------------------------------------------------------------
203  CALL fine2coarse(xdata,ydata,x,y,out,msk=NINT(glace01)==1)
204  frac_ice=REAL(out)
205
206END SUBROUTINE sea_ice
207!
208!-------------------------------------------------------------------------------
209
210
211!-------------------------------------------------------------------------------
212!
213SUBROUTINE rugsoro(xrel, yrel, relief, xmod, ymod, rugs)
214!
215!-------------------------------------------------------------------------------
216! Author:  Z.X. Li (april 1st 1994) ; Modifications: august 23rd 1995.
217! Purpose: Compute rugosity due to orography by using standard dev in a 1x1 cell
218!-------------------------------------------------------------------------------
219  IMPLICIT NONE
220!-------------------------------------------------------------------------------
221! Arguments:
222  REAL, INTENT(IN)  :: xrel(:), yrel(:)        !--- INPUT  FIELD X AND Y COORD.
223  REAL, INTENT(IN)  :: relief(:,:)             !--- INPUT  FIELD
224  REAL, INTENT(IN)  :: xmod(:), ymod(:)        !--- OUTPUT FIELD X AND Y COORD.
225  REAL, INTENT(OUT) :: rugs(SIZE(xmod),SIZE(ymod)) !--- OUTPUT FIELD
226!-------------------------------------------------------------------------------
227! Local variable:
228  INTEGER           :: k, nn
229  INTEGER, PARAMETER:: itmp=360, jtmp=180
230  DOUBLE PRECISION  :: out(SIZE(xmod),SIZE(xmod)), amin, amax
231  DOUBLE PRECISION  :: cham1tmp(itmp,jtmp), cham2tmp(itmp,jtmp)
232  REAL              :: xtmp(itmp), ytmp(jtmp)
233!-------------------------------------------------------------------------------
234
235!--- TEST INPUT FILE FITS FOR ONE DEGREE RESOLUTION
236  xtmp(:)=4.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(itmp),k=1,itmp)]
237  ytmp(:)=2.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(jtmp),k=1,jtmp)]
238  CALL fine2coarse(xrel,yrel,xtmp,ytmp,cham1tmp,relief,d_o2=cham2tmp)
239  cham2tmp(:,:)=cham2tmp(:,:)-cham1tmp(:,:)**2
240  nn=COUNT(cham2tmp<=-7.5)
241  IF(nn/=0) THEN
242    PRINT*,'Problem for rugsoro ; std**2 < -7.5 for several points: ',nn
243    CALL ABORT_GCM("", "", 1)
244  END IF
245  nn=COUNT(cham2tmp<0)
246  IF(nn/=0) PRINT*,'Problem for rugsoro ; std**2 < 0. for several points: ',nn
247  WHERE(cham2tmp<0.0) cham2tmp=0.0d0
248  cham2tmp(:,:)=SQRT(cham2tmp(:,:))
249  amin=MINVAL(cham2tmp); amax=MAXVAL(cham2tmp)
250  PRINT*, 'Ecart-type 1x1:', amin, amax
251
252!--- COMPUTE RUGOSITY AT REQUIRED SCALE
253  WHERE(cham2tmp<0.001d0) cham2tmp=0.001d0
254  CALL fine2coarse(xtmp,ytmp,xmod,ymod,out,REAL(LOG(cham2tmp)))
255  out=EXP(out)
256  amin=MINVAL(out); amax=MAXVAL(out)
257  PRINT*, 'Ecart-type du modele:', amin, amax
258  out=out/amax*20.0d0
259  amin=MINVAL(out); amax=MAXVAL(out)
260  PRINT*, 'Longueur de rugosite du modele:', amin, amax
261  rugs=REAL(out)
262
263END SUBROUTINE rugsoro
264!
265!-------------------------------------------------------------------------------
266
267
268!-------------------------------------------------------------------------------
269!
270SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,distance)
271!
272!-------------------------------------------------------------------------------
273! Author:  Laurent Li (december 30th 1996).
274! Purpose: Compute min. distance (along big circle) between 2 points in degrees.
275!-------------------------------------------------------------------------------
276  IMPLICIT NONE
277!-------------------------------------------------------------------------------
278! Arguments:
279  REAL, INTENT(IN) :: rf_lon, rf_lat  !--- Reference point coordinates (degrees)
280  REAL, INTENT(IN) :: rlon(:), rlat(:)!--- Points longitudes/latitudes (degrees)
281  REAL, INTENT(OUT):: distance(SIZE(rlon),SIZE(rlat)) !--- Distance    (degrees)
282!-------------------------------------------------------------------------------
283! Local variables:
284  LOGICAL, SAVE :: first=.TRUE.
285  REAL    :: pa, pb, cpa, cpab, spa, spab, crlo(SIZE(rlon))
286  INTEGER :: i, j
287!-------------------------------------------------------------------------------
288  IF(first) THEN
289    pi=4.0*ATAN(1.0); deg2rad=pi/180.0; first=.FALSE.
290  END IF
291  crlo(:)=COS((rf_lon-rlon(:))*deg2rad)     !--- COS of points 1 and 2 angle
292  pa=(90.0-rf_lat)*deg2rad                  !--- North Pole - Point 1 distance
293  cpa=COS(pa); spa=SIN(pa)
294  DO j=1,SIZE(rlat)
295    pb=(90.0-rlat(j))*deg2rad               !--- North Pole - Point 2 distance
296    cpab=cpa*COS(pb); spab=spa*SIN(pb)
297    distance(:,j)=ACOS(cpab+spab*crlo(:))/deg2rad
298  END DO
299
300END SUBROUTINE dist_sphe
301!
302!-------------------------------------------------------------------------------
303
304END MODULE grid_atob_m
305!
306!*******************************************************************************
307
Note: See TracBrowser for help on using the repository browser.