1 | PROGRAM create_coordinates |
---|
2 | !!----------------------------------------------------------- |
---|
3 | !! |
---|
4 | !! to create regional coordinates file (e.g. 1_coordinates.nc) |
---|
5 | !! for fine grid from coarse grid (e.g. coordinates.nc) |
---|
6 | !! |
---|
7 | !! to make it we use a 4th order polynomial interpolation |
---|
8 | !! |
---|
9 | !! Created by Brice Lemaire on 12/2009. |
---|
10 | !! |
---|
11 | !!----------------------------------------------------------- |
---|
12 | USE netcdf |
---|
13 | USE cfg_tools |
---|
14 | USE mixed_grid |
---|
15 | USE domain |
---|
16 | !------------ |
---|
17 | IMPLICIT NONE |
---|
18 | ! |
---|
19 | INTEGER :: narg,iargc |
---|
20 | INTEGER :: nstat, imodulo |
---|
21 | INTEGER :: nik, njk, npk |
---|
22 | INTEGER :: ival |
---|
23 | CHARACTER(len=80) :: cnmlname, cfingrdname |
---|
24 | ! |
---|
25 | ! |
---|
26 | ! |
---|
27 | !************************************* |
---|
28 | ! Read input file (namelist file) |
---|
29 | ! using types.f90 |
---|
30 | !************************************* |
---|
31 | narg = iargc() |
---|
32 | ! |
---|
33 | IF (narg == 0) THEN |
---|
34 | cnmlname = 'namelist.input' |
---|
35 | ELSE |
---|
36 | CALL getarg(1,cnmlname) |
---|
37 | ENDIF |
---|
38 | ! |
---|
39 | CALL read_namelist(cnmlname) |
---|
40 | ! |
---|
41 | nstat = read_coordinates(TRIM(cn_parent_coordinate_file),scoagrd) |
---|
42 | IF (nstat/=1) THEN |
---|
43 | WRITE(*,*) 'unable to open netcdf file : ',cn_parent_coordinate_file |
---|
44 | STOP |
---|
45 | ENDIF |
---|
46 | ! |
---|
47 | nsizex = SIZE(scoagrd%glamt,1) |
---|
48 | nsizey = SIZE(scoagrd%glamt,2) |
---|
49 | ! |
---|
50 | WRITE(*,*) '' |
---|
51 | WRITE(*,*) 'Size of input matrix: ' |
---|
52 | WRITE(*,*) '(',nsizex,';', nsizey,')' |
---|
53 | WRITE(*,*) '' |
---|
54 | ! |
---|
55 | WRITE(*,*) 'Domain: ' |
---|
56 | WRITE(*,*) '' |
---|
57 | WRITE(*,*) ' ', ' min(1,1:nsizey) ', ' max(nsizex,1:nsizey) ' |
---|
58 | WRITE(*,*) 'longitude: ', MINVAL(scoagrd%glamt(1,1:nsizey)), ' --> ', MAXVAL(scoagrd%glamf(nsizex,1:nsizey)) |
---|
59 | WRITE(*,*) 'latitude: ', MINVAL(scoagrd%gphit(1:nsizex,1)), ' --> ', MAXVAL(scoagrd%gphif(1:nsizex,nsizey)) |
---|
60 | WRITE(*,*) '' |
---|
61 | ! |
---|
62 | !************************************* |
---|
63 | ! Check the values read in the namelist |
---|
64 | !************************************* |
---|
65 | IF(nn_imin.LE.0.OR.nn_imin.GT.nsizex)THEN |
---|
66 | WRITE(*,*) 'Wrong values in namelist:' |
---|
67 | WRITE(*,*) 'nn_imin must be greater than 0 and less than', nsizex |
---|
68 | STOP |
---|
69 | ENDIF |
---|
70 | ! |
---|
71 | IF(nn_imax.LE.0.OR.nn_imax.GT.nsizex)THEN |
---|
72 | WRITE(*,*) 'Wrong values in namelist:' |
---|
73 | WRITE(*,*) 'nn_imax must be greater than 0 and less than', nsizex |
---|
74 | STOP |
---|
75 | ENDIF |
---|
76 | ! |
---|
77 | IF(nn_jmin.LE.0.OR.nn_jmin.GT.nsizey)THEN |
---|
78 | WRITE(*,*) 'Wrong values in namelist:' |
---|
79 | WRITE(*,*) 'nn_jmin must be greater than 0 and less than', nsizey |
---|
80 | STOP |
---|
81 | ENDIF |
---|
82 | ! |
---|
83 | IF(nn_jmax.LE.0.OR.nn_jmax.GT.nsizey)THEN |
---|
84 | WRITE(*,*) 'Wrong values in namelist:' |
---|
85 | WRITE(*,*) 'nn_jmax must be greater than 0 and less than', nsizey |
---|
86 | STOP |
---|
87 | ENDIF |
---|
88 | ! |
---|
89 | !************************************* |
---|
90 | ! Define the sub-domain chosen by user |
---|
91 | !************************************* |
---|
92 | WRITE(*,*) 'Domain defined by user: ' |
---|
93 | WRITE(*,*) '' |
---|
94 | WRITE(*,*) ' ', ' min ', ' max ' |
---|
95 | WRITE(*,*) 'longitude: ', MINVAL(scoagrd%glamt(nn_imin,nn_jmin:nn_jmax)), ' --> ', MAXVAL(scoagrd%glamf(nn_imax,nn_jmin:nn_jmax)) |
---|
96 | WRITE(*,*) 'latitude: ', MINVAL(scoagrd%gphit(nn_imin:nn_imax,nn_jmin)), ' --> ', MAXVAL(scoagrd%gphif(nn_imin:nn_imax,nn_jmax)) |
---|
97 | WRITE(*,*) '' |
---|
98 | ! |
---|
99 | IF(nn_imin.EQ.nn_imax.AND.nn_jmin.EQ.nn_jmax) THEN |
---|
100 | nglobal = .TRUE. |
---|
101 | WRITE(*,*) 'Size of domain: GLOBAL' |
---|
102 | ELSE |
---|
103 | nglobal = .FALSE. |
---|
104 | WRITE(*,*) 'Size of domain: REGIONAL' |
---|
105 | ENDIF |
---|
106 | ! |
---|
107 | IF(cn_position_pivot.EQ.'F-grid') THEN !case of ORCA05 & ORCA025 grids |
---|
108 | npivot = 0 |
---|
109 | ELSEIF(cn_position_pivot.EQ.'T-grid') THEN !case of ORCA2 grid |
---|
110 | npivot = 1 |
---|
111 | ENDIF |
---|
112 | ! |
---|
113 | nmid = nsizex/2 + npivot |
---|
114 | ! |
---|
115 | ! |
---|
116 | ! |
---|
117 | !************************************* |
---|
118 | ! Redefine for particular cases |
---|
119 | !************************************* |
---|
120 | IF(((nn_imax - nn_imin).EQ.-1).OR.((nn_imax - nn_imin).EQ.-2)) THEN |
---|
121 | nn_imax = nn_imin |
---|
122 | ELSEIF(ABS(nn_imin - nn_imax).GE.nsizex-1) THEN |
---|
123 | nn_imax = nn_imin |
---|
124 | ELSEIF((nn_imin.EQ.1).AND.(nn_imax.EQ.nsizex)) THEN |
---|
125 | nn_imax = nn_imin |
---|
126 | ELSEIF(nn_imin.EQ.1.AND.nn_imax.GT.nn_imin) THEN |
---|
127 | nn_imin = nsizey-2 |
---|
128 | ELSEIF(nn_imax.EQ.nsizex) THEN |
---|
129 | nn_imax = 3 |
---|
130 | ELSEIF(nn_imin.EQ.nsizex) THEN |
---|
131 | nn_imin = nsizex-1 |
---|
132 | ELSEIF(nn_imax.EQ.1.AND.nn_imax.NE.nn_imin) THEN |
---|
133 | nn_imax = 2 |
---|
134 | ENDIF |
---|
135 | ! |
---|
136 | ! |
---|
137 | ! |
---|
138 | !************************************* |
---|
139 | ! We want to fix equator along T and U-points |
---|
140 | !************************************* |
---|
141 | imodulo = MOD(nn_rhoy,2) |
---|
142 | ! |
---|
143 | IF(.NOT.nglobal.AND.(scoagrd%gphit(nn_imin,nn_jmin)*scoagrd%gphit(nn_imin,nn_jmax)).LT.0.AND.imodulo.EQ.0) THEN |
---|
144 | nequator = 1 |
---|
145 | ELSEIF(nglobal.AND.imodulo.EQ.0) THEN |
---|
146 | nequator = 1 |
---|
147 | ELSE |
---|
148 | nequator = 0 |
---|
149 | ENDIF |
---|
150 | ! |
---|
151 | ! |
---|
152 | ! |
---|
153 | !************************************* |
---|
154 | ! CREATE MIXED GRID |
---|
155 | !************************************* |
---|
156 | CALL define_domain |
---|
157 | ! |
---|
158 | CALL define_mixed_grid |
---|
159 | ! |
---|
160 | ! |
---|
161 | ! |
---|
162 | !************************************* |
---|
163 | !!! CALCULATE FINE GRID DIMENSIONS |
---|
164 | !************************************* |
---|
165 | IF(.NOT.nglobal) THEN |
---|
166 | IF(nn_rhox.EQ.1) THEN |
---|
167 | nxfine = nxcoag |
---|
168 | ELSE |
---|
169 | nxfine = (nxcoag-2)*nn_rhox + 1 |
---|
170 | ENDIF |
---|
171 | ! |
---|
172 | IF(nn_rhoy.EQ.1) THEN |
---|
173 | nyfine = nycoag |
---|
174 | ELSE |
---|
175 | nyfine = (nycoag-2)*nn_rhoy + 1 |
---|
176 | ENDIF |
---|
177 | ! |
---|
178 | ELSEIF(nglobal) THEN |
---|
179 | IF(nn_rhox.GT.1) THEN |
---|
180 | nxfine = (nxgmix - 4*(nn_rhox-1))/2 |
---|
181 | ELSEIF(nn_rhox.EQ.1) THEN |
---|
182 | nxfine = nsizex |
---|
183 | ENDIF |
---|
184 | ! |
---|
185 | IF(nn_rhoy.GT.1) THEN |
---|
186 | nyfine = (nygmix - (2*nn_rhoy))/2 - nn_rhoy + nequator |
---|
187 | ELSEIF(nn_rhoy.EQ.1) THEN |
---|
188 | nyfine = nsizey |
---|
189 | ENDIF |
---|
190 | ! |
---|
191 | ENDIF |
---|
192 | ! |
---|
193 | WRITE(*,*) '' |
---|
194 | WRITE(*,*) '*** SIZE OF FINE GRID ***' |
---|
195 | WRITE(*,*) nxfine, ' x ', nyfine |
---|
196 | WRITE(*,*) '' |
---|
197 | ! |
---|
198 | ! |
---|
199 | ! |
---|
200 | !************************************* |
---|
201 | ! Interpolation inside the mixed grid |
---|
202 | ! from cfg_tools.f90 |
---|
203 | !************************************* |
---|
204 | IF(nn_rhox.GT.1.OR.nn_rhoy.GT.1) THEN |
---|
205 | CALL interp_grid |
---|
206 | ENDIF |
---|
207 | ! |
---|
208 | ! |
---|
209 | ! |
---|
210 | !************************************* |
---|
211 | ! Define name of child coordinate file |
---|
212 | ! coordinates.nc -> 1_coordinates.nc |
---|
213 | !************************************* |
---|
214 | CALL set_child_name(cn_parent_coordinate_file,cfingrdname) |
---|
215 | ! |
---|
216 | ! |
---|
217 | ! |
---|
218 | !************************************* |
---|
219 | ! Allocation of child grid elements |
---|
220 | ! from types.f90 |
---|
221 | !************************************* |
---|
222 | CALL grid_allocate(sfingrd,nxfine,nyfine) |
---|
223 | ! |
---|
224 | ! |
---|
225 | ! |
---|
226 | !************************************* |
---|
227 | ! Break the mixed grid smixgrd into 4 grids |
---|
228 | ! from cfg_tools.f90 |
---|
229 | !************************************* |
---|
230 | CALL child_grid |
---|
231 | ! |
---|
232 | ! |
---|
233 | ! |
---|
234 | !************************************* |
---|
235 | ! Read parent coordinate file |
---|
236 | ! from readwrite.f90 |
---|
237 | !************************************* |
---|
238 | nstat = write_coordinates(cfingrdname,sfingrd,nxfine,nyfine) |
---|
239 | IF (nstat/=1) THEN |
---|
240 | WRITE(*,*)"unable to write netcdf file : ",cfingrdname |
---|
241 | STOP |
---|
242 | ENDIF |
---|
243 | ! |
---|
244 | CALL grid_deallocate(scoagrd) |
---|
245 | CALL grid_deallocate(sfingrd) |
---|
246 | CALL mixed_grid_deallocate(smixgrd) |
---|
247 | ! |
---|
248 | ! |
---|
249 | ! |
---|
250 | END PROGRAM create_coordinates |
---|