source: utils/tools/NESTING/src/agrif_create_coordinates.f90 @ 9632

Last change on this file since 9632 was 9632, checked in by jchanut, 3 years ago

Correct child velocity points lat/lon interpolation, #2082

  • Property svn:keywords set to Id
File size: 6.2 KB
Line 
1PROGRAM create_coordinate
2  !
3  USE NETCDF
4  USE agrif_readwrite
5  USE agrif_interpolation
6  !       
7  !****************************************************************     
8  !                        *
9  !************************************************************************
10  ! Fortran 95 OPA Nesting tools                *
11  !                           *
12  !     Copyright (C) 2005 Florian Lemarié (Florian.Lemarie@imag.fr) *
13  !                           *
14  !************************************************************************
15  !
16  ! PROGRAM CREATE_COORDINATE                *
17  !                        *
18  ! program to create coordinates file for child grid    *
19  ! (child grid defined by imin,imax,jmin,jmax,rho)      *
20  !                        *
21  !****************************************************************
22  !
23  ! variables declaration
24  !
25  REAL*8 :: rpi,ra,rad
26  CHARACTER*100 :: Child_filename 
27  INTEGER :: i
28  TYPE(Coordinates) :: G0,G1
29  INTEGER :: narg,iargc
30  CHARACTER(len=80) :: namelistname
31
32  narg = iargc()
33
34  IF (narg == 0) THEN
35     namelistname = 'namelist.input'
36  ELSE
37     CALL getarg(1,namelistname)
38  ENDIF
39  !
40  ! read input file (namelist.input)
41  !
42  CALL read_namelist(namelistname)
43  !
44
45  !     
46  ! read parent coodinates file
47  !
48  status = Read_Coordinates(TRIM(parent_coordinate_file),G0)
49  !
50  ! define name of child coordinate file
51  !     
52  CALL set_child_name(parent_coordinate_file,Child_filename) 
53  !
54  IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN                   
55     WRITE(*,*) 'ERROR ***** bad child grid definition ...'
56     WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
57     WRITE(*,*) ' '
58     STOP
59  ENDIF
60  !                     
61  WRITE(*,*) ' '
62  WRITE(*,*)'Size of the High resolution grid: ',nxfin,' x ',nyfin
63  WRITE(*,*) ' '
64  !
65  ! allocation of child grid elements
66  !
67  CALL agrif_grid_allocate(G1,nxfin,nyfin)
68  !
69  !
70  !  check potential longitude problems
71  !     
72  IF( G0%glamt(imin,jmin) > G0%glamt(imax,jmax) ) THEN           
73     WRITE(*,*) '    '
74     WHERE ( G0%glamt < 0 )
75        G0%glamt = G0%glamt + 360.
76     END WHERE
77     WHERE ( G0%glamf < 0 )
78        G0%glamf = G0%glamf + 360.
79     END WHERE
80  ENDIF
81
82  !     
83  ! interpolation from parent grid to child grid for
84  ! T points (cell center)
85  ! F points (lower left corner)
86  ! U points (cell left side)
87  ! V points (cell top side)
88  !    glam = longitude
89  !    gphi = latitude
90  !
91
92  !     
93  !> M. Dunphy ticket 2082:
94  CALL agrif_interp(G0%glamt,G1%glamt,'T')
95  CALL agrif_interp(G0%glamf,G1%glamf,'F')
96!  G1%glamu = G1%glamf
97!  G1%glamv = G1%glamt
98  CALL agrif_interp(G0%glamu,G1%glamu,'U')
99  CALL agrif_interp(G0%glamv,G1%glamv,'V') 
100  !
101  CALL agrif_interp(G0%gphit,G1%gphit,'T')
102  CALL agrif_interp(G0%gphif,G1%gphif,'F')
103!  G1%gphiu = G1%gphit
104!  G1%gphiv = G1%gphif
105  CALL agrif_interp(G0%gphiu,G1%gphiu,'U')
106  CALL agrif_interp(G0%gphiv,G1%gphiv,'V')
107  !<  M. Dunphy ticket 2082
108  !
109  !
110  rpi = 4.*ATAN(1.)
111  rad = rpi/180.
112  ra  = 6371229.
113  !
114  ! Compute scale factors e1 e2
115  !
116  DO j=1,nyfin
117     DO i=2,nxfin
118        G1%e1t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamu(i,j)-G1%glamu(i-1,j)))**2 &
119             + (G1%gphiu(i,j)-G1%gphiu(i-1,j))**2)
120        G1%e1v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamf(i,j)-G1%glamf(i-1,j)))**2 &
121             + (G1%gphif(i,j)-G1%gphif(i-1,j))**2)                                                             
122     END DO
123  END DO
124  !     
125  G1%e1t(1,:)=G1%e1t(2,:)
126  G1%e1v(1,:)=G1%e1v(2,:)
127  !
128  DO j=1,nyfin
129     DO i=1,nxfin-1
130        G1%e1u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamt(i+1,j)-G1%glamt(i,j)))**2 &
131             + (G1%gphit(i+1,j)-G1%gphit(i,j))**2)
132        G1%e1f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamv(i+1,j)-G1%glamv(i,j)))**2 &
133             + (G1%gphiv(i+1,j)-G1%gphiv(i,j))**2)                                                           
134     END DO
135  END DO
136  !
137  G1%e1u(nxfin,:)=G1%e1u(nxfin-1,:)
138  G1%e1f(nxfin,:)=G1%e1f(nxfin-1,:)
139  !
140  DO j=2,nyfin
141     DO i=1,nxfin                                     
142        G1%e2t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamv(i,j)-G1%glamv(i,j-1)))**2 &
143             + (G1%gphiv(i,j)-G1%gphiv(i,j-1))**2)
144        G1%e2u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamf(i,j)-G1%glamf(i,j-1)))**2 &
145             + (G1%gphif(i,j)-G1%gphif(i,j-1))**2)                                                               
146     END DO
147  END DO
148  !
149  G1%e2t(:,1)=G1%e2t(:,2)
150  G1%e2u(:,1)=G1%e2u(:,2)
151  !
152  DO j=1,nyfin-1
153     DO i=1,nxfin
154        G1%e2v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamt(i,j+1)-G1%glamt(i,j)))**2 &
155             + (G1%gphit(i,j+1)-G1%gphit(i,j))**2)
156        G1%e2f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamu(i,j+1)-G1%glamu(i,j)))**2 &
157             + (G1%gphiu(i,j+1)-G1%gphiu(i,j))**2)
158     END DO
159  END DO
160  !
161  G1%e2v(:,nyfin)=G1%e2v(:,nyfin-1)
162  G1%e2f(:,nyfin)=G1%e2f(:,nyfin-1)
163
164
165  CALL agrif_interp(G0%e1t,G1%e1t,'T')
166  G1%e1t = G1%e1t / REAL(irafx)
167  CALL agrif_interp(G0%e2t,G1%e2t,'T')
168  G1%e2t = G1%e2t / REAL(irafy)
169
170  CALL agrif_interp(G0%e1u,G1%e1u,'U')
171  G1%e1u = G1%e1u / REAL(irafx)
172  CALL agrif_interp(G0%e2u,G1%e2u,'U')
173  G1%e2u = G1%e2u / REAL(irafy) 
174
175  CALL agrif_interp(G0%e1v,G1%e1v,'V')
176  G1%e1v = G1%e1v / REAL(irafx)
177  CALL agrif_interp(G0%e2v,G1%e2v,'V')
178  G1%e2v = G1%e2v / REAL(irafy) 
179
180  CALL agrif_interp(G0%e1f,G1%e1f,'F')
181  G1%e1f = G1%e1f / REAL(irafx)
182  CALL agrif_interp(G0%e2f,G1%e2f,'F')
183  G1%e2f = G1%e2f / REAL(irafy)               
184
185  !
186  WHERE ( G1%glamt > 180 )
187     G1%glamt = G1%glamt - 360.
188  END WHERE
189  WHERE ( G1%glamf > 180 )
190     G1%glamf = G1%glamf - 360.
191  END WHERE
192  WHERE ( G1%glamu > 180 )
193     G1%glamu = G1%glamu - 360.
194  END WHERE
195  WHERE ( G1%glamv > 180 )
196     G1%glamv = G1%glamv - 360.
197  END WHERE
198  !
199  !
200  G1%nav_lon=G1%glamt
201  G1%nav_lat=G1%gphit 
202  !             
203  !
204  ! Write interpolation result in child coodinates file
205  !     
206  status = Write_Coordinates(TRIM(Child_filename),G1)
207
208  !
209  WRITE(*,*) 'Child domain position : '
210  WRITE(*,*) 'latmin =',G1%gphit(3,3)
211  WRITE(*,*) 'latmax =',G1%gphit(nxfin-2,nyfin-2)
212  WRITE(*,*) 'lonmin =',G1%glamt(3,3)
213  WRITE(*,*) 'lonmax =',G1%glamt(nxfin-2,nyfin-2)
214  STOP
215END PROGRAM create_coordinate
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
Note: See TracBrowser for help on using the repository browser.