New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_create_coordinates.f90 in branches/UKMO/restart_datestamp/NEMOGCM/TOOLS/NESTING/src – NEMO

source: branches/UKMO/restart_datestamp/NEMOGCM/TOOLS/NESTING/src/agrif_create_coordinates.f90 @ 5420

Last change on this file since 5420 was 5420, checked in by davestorkey, 9 years ago

Remove keyword updating from UKMO restart_datestamp branch.

File size: 6.0 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  CALL agrif_interp(G0%glamt,G1%glamt,'T')           
94  CALL agrif_interp(G0%glamf,G1%glamf,'F')       
95  G1%glamu = G1%glamf
96  G1%glamv = G1%glamt   
97  !
98  CALL agrif_interp(G0%gphit,G1%gphit,'T')
99  CALL agrif_interp(G0%gphif,G1%gphif,'F')
100  G1%gphiu = G1%gphit
101  G1%gphiv = G1%gphif
102  !
103  !
104  rpi = 4.*ATAN(1.)
105  rad = rpi/180.
106  ra  = 6371229.
107  !
108  ! Compute scale factors e1 e2
109  !
110  DO j=1,nyfin
111     DO i=2,nxfin
112        G1%e1t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamu(i,j)-G1%glamu(i-1,j)))**2 &
113             + (G1%gphiu(i,j)-G1%gphiu(i-1,j))**2)
114        G1%e1v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamf(i,j)-G1%glamf(i-1,j)))**2 &
115             + (G1%gphif(i,j)-G1%gphif(i-1,j))**2)                                                             
116     END DO
117  END DO
118  !     
119  G1%e1t(1,:)=G1%e1t(2,:)
120  G1%e1v(1,:)=G1%e1v(2,:)
121  !
122  DO j=1,nyfin
123     DO i=1,nxfin-1
124        G1%e1u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamt(i+1,j)-G1%glamt(i,j)))**2 &
125             + (G1%gphit(i+1,j)-G1%gphit(i,j))**2)
126        G1%e1f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamv(i+1,j)-G1%glamv(i,j)))**2 &
127             + (G1%gphiv(i+1,j)-G1%gphiv(i,j))**2)                                                           
128     END DO
129  END DO
130  !
131  G1%e1u(nxfin,:)=G1%e1u(nxfin-1,:)
132  G1%e1f(nxfin,:)=G1%e1f(nxfin-1,:)
133  !
134  DO j=2,nyfin
135     DO i=1,nxfin                                     
136        G1%e2t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamv(i,j)-G1%glamv(i,j-1)))**2 &
137             + (G1%gphiv(i,j)-G1%gphiv(i,j-1))**2)
138        G1%e2u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamf(i,j)-G1%glamf(i,j-1)))**2 &
139             + (G1%gphif(i,j)-G1%gphif(i,j-1))**2)                                                               
140     END DO
141  END DO
142  !
143  G1%e2t(:,1)=G1%e2t(:,2)
144  G1%e2u(:,1)=G1%e2u(:,2)
145  !
146  DO j=1,nyfin-1
147     DO i=1,nxfin
148        G1%e2v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamt(i,j+1)-G1%glamt(i,j)))**2 &
149             + (G1%gphit(i,j+1)-G1%gphit(i,j))**2)
150        G1%e2f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamu(i,j+1)-G1%glamu(i,j)))**2 &
151             + (G1%gphiu(i,j+1)-G1%gphiu(i,j))**2)
152     END DO
153  END DO
154  !
155  G1%e2v(:,nyfin)=G1%e2v(:,nyfin-1)
156  G1%e2f(:,nyfin)=G1%e2f(:,nyfin-1)
157
158
159  CALL agrif_interp(G0%e1t,G1%e1t,'T')
160  G1%e1t = G1%e1t / REAL(irafx)
161  CALL agrif_interp(G0%e2t,G1%e2t,'T')
162  G1%e2t = G1%e2t / REAL(irafy)
163
164  CALL agrif_interp(G0%e1u,G1%e1u,'U')
165  G1%e1u = G1%e1u / REAL(irafx)
166  CALL agrif_interp(G0%e2u,G1%e2u,'U')
167  G1%e2u = G1%e2u / REAL(irafy) 
168
169  CALL agrif_interp(G0%e1v,G1%e1v,'V')
170  G1%e1v = G1%e1v / REAL(irafx)
171  CALL agrif_interp(G0%e2v,G1%e2v,'V')
172  G1%e2v = G1%e2v / REAL(irafy) 
173
174  CALL agrif_interp(G0%e1f,G1%e1f,'F')
175  G1%e1f = G1%e1f / REAL(irafx)
176  CALL agrif_interp(G0%e2f,G1%e2f,'F')
177  G1%e2f = G1%e2f / REAL(irafy)               
178
179  !
180  WHERE ( G1%glamt > 180 )
181     G1%glamt = G1%glamt - 360.
182  END WHERE
183  WHERE ( G1%glamf > 180 )
184     G1%glamf = G1%glamf - 360.
185  END WHERE
186  WHERE ( G1%glamu > 180 )
187     G1%glamu = G1%glamu - 360.
188  END WHERE
189  WHERE ( G1%glamv > 180 )
190     G1%glamv = G1%glamv - 360.
191  END WHERE
192  !
193  !
194  G1%nav_lon=G1%glamt
195  G1%nav_lat=G1%gphit 
196  !             
197  !
198  ! Write interpolation result in child coodinates file
199  !     
200  status = Write_Coordinates(TRIM(Child_filename),G1)
201
202  !
203  WRITE(*,*) 'Child domain position : '
204  WRITE(*,*) 'latmin =',G1%gphit(3,3)
205  WRITE(*,*) 'latmax =',G1%gphit(nxfin-2,nyfin-2)
206  WRITE(*,*) 'lonmin =',G1%glamt(3,3)
207  WRITE(*,*) 'lonmax =',G1%glamt(nxfin-2,nyfin-2)
208  STOP
209END PROGRAM create_coordinate
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
Note: See TracBrowser for help on using the repository browser.