source: LMDZ5/trunk/libf/dynlonlat_phylonlat/grid_atob_m.f90 @ 2360

Last change on this file since 2360 was 2299, checked in by dcugnet, 9 years ago

In dyn3d/:
etat0dyn_netcdf.F90: "startget_dyn3d" syntax slightly simplified.
dynredem.F90: Shortcut routines (put_var*, cre_var,
dynredem_write_*, dynredem_read_u)

modified to match dyn3dmem version and put in

module dyredem_mod.F90.
dynetat0.F90 -> *.f90: Few simplifications (no usage of NC_DOUBLE
needed => no precompilation)

Add tracers initialization in the isotope case

suppressed by accident.
dynredem_mod.F90: Created to mimic dyn3dmem equivalent.

In dyn3dmem/:
dynetat0_loc.F -> *.f90: Converted into fortran 90 to match the dyn3d
version.
dynredem_loc.F -> *.F90: Converted into fortran 90.
dynredem_mod.F90: Add some shortcut routines to match the dyn3d
version.

In phylmd/:
phyredem.F90: Bug fix: nsw instead of nsoilmx was used as
Tsoil second maximum index.

Bug fix: fevap instead of snow was saved for

"SNOW".
etat0phys_netcdf.F90: "filtreg_mod" module usage suppressed.

Local variable rugo computation removed (not

used).

In dynlonlat_phylonlat/:
grid_atob_m.F90 -> *.f90 DOUBLE PRECISION variables usage removed.

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