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

Last change on this file since 10025 was 10025, checked in by clem, 2 years ago

nesting tools are partly rewritten (mostly for create_coordinates and bathy) to get better functionality. Now you can use the nesting to either define an agrif zoom or a regional domain (for bdy purposes). Also, the nesting tools output a domain_cfg.nc that can be directly used in NEMO4 without the need of DOMAINcfg tool. The option of median average for bathymetry interpolation still does not work properly but it's not new

  • Property svn:keywords set to Id
File size: 6.1 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  CALL read_namelist(namelistname)
42  !     
43  ! read parent coodinates file
44  status = Read_Coordinates(TRIM(parent_coordinate_file),G0)
45  !
46  ! define name of child coordinate file
47  CALL set_child_name(parent_coordinate_file,Child_filename) 
48  !
49  IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN                   
50     WRITE(*,*) 'ERROR ***** bad child grid definition ...'
51     WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
52     WRITE(*,*) ' '
53     STOP
54  ENDIF
55  !                     
56  WRITE(*,*) ' '
57  WRITE(*,*)'Size of the High resolution grid: ',nxfin,' x ',nyfin
58  WRITE(*,*) ' '
59  !
60  ! allocation of child grid elements
61  CALL agrif_grid_allocate(G1,nxfin,nyfin)
62  !
63  !  check potential longitude problems
64  IF( G0%glamt(imin,jmin) > G0%glamt(imax,jmax) ) THEN           
65     WHERE ( G0%glamt < 0 )   G0%glamt = G0%glamt + 360.
66     WHERE ( G0%glamf < 0 )   G0%glamf = G0%glamf + 360.
67  ENDIF
68  !     
69  ! interpolation from parent grid to child grid for
70  ! T points (cell center)
71  ! F points (lower left corner)
72  ! U points (cell left side)
73  ! V points (cell top side)
74  !    glam = longitude
75  !    gphi = latitude
76  !
77  !     
78  CALL agrif_interp(G0%glamt,G1%glamt,'T')
79  CALL agrif_interp(G0%glamf,G1%glamf,'F')
80  CALL agrif_interp(G0%glamu,G1%glamu,'U')
81  CALL agrif_interp(G0%glamv,G1%glamv,'V') 
82  !
83  CALL agrif_interp(G0%gphit,G1%gphit,'T')
84  CALL agrif_interp(G0%gphif,G1%gphif,'F')
85  CALL agrif_interp(G0%gphiu,G1%gphiu,'U')
86  CALL agrif_interp(G0%gphiv,G1%gphiv,'V')
87  !
88  !
89  rpi = 4.*ATAN(1.)
90  rad = rpi/180.
91  ra  = 6371229.
92  !
93  ! Compute scale factors e1 e2
94!  DO j=1,nyfin
95!     DO i=2,nxfin
96!        G1%e1t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamu(i,j)-G1%glamu(i-1,j)))**2 &
97!             + (G1%gphiu(i,j)-G1%gphiu(i-1,j))**2)
98!        G1%e1v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamf(i,j)-G1%glamf(i-1,j)))**2 &
99!             + (G1%gphif(i,j)-G1%gphif(i-1,j))**2)                                                             
100!     END DO
101!  END DO
102!  !     
103!  G1%e1t(1,:)=G1%e1t(2,:)
104!  G1%e1v(1,:)=G1%e1v(2,:)
105!  !
106!  DO j=1,nyfin
107!     DO i=1,nxfin-1
108!        G1%e1u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamt(i+1,j)-G1%glamt(i,j)))**2 &
109!             + (G1%gphit(i+1,j)-G1%gphit(i,j))**2)
110!        G1%e1f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamv(i+1,j)-G1%glamv(i,j)))**2 &
111!             + (G1%gphiv(i+1,j)-G1%gphiv(i,j))**2)                                                           
112!     END DO
113!  END DO
114!  !
115!  G1%e1u(nxfin,:)=G1%e1u(nxfin-1,:)
116!  G1%e1f(nxfin,:)=G1%e1f(nxfin-1,:)
117!  !
118!  DO j=2,nyfin
119!     DO i=1,nxfin                                     
120!        G1%e2t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamv(i,j)-G1%glamv(i,j-1)))**2 &
121!             + (G1%gphiv(i,j)-G1%gphiv(i,j-1))**2)
122!        G1%e2u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamf(i,j)-G1%glamf(i,j-1)))**2 &
123!             + (G1%gphif(i,j)-G1%gphif(i,j-1))**2)                                                               
124!     END DO
125!  END DO
126!  !
127!  G1%e2t(:,1)=G1%e2t(:,2)
128!  G1%e2u(:,1)=G1%e2u(:,2)
129!  !
130!  DO j=1,nyfin-1
131!     DO i=1,nxfin
132!        G1%e2v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamt(i,j+1)-G1%glamt(i,j)))**2 &
133!             + (G1%gphit(i,j+1)-G1%gphit(i,j))**2)
134!        G1%e2f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamu(i,j+1)-G1%glamu(i,j)))**2 &
135!             + (G1%gphiu(i,j+1)-G1%gphiu(i,j))**2)
136!     END DO
137!  END DO
138!  !
139!  G1%e2v(:,nyfin)=G1%e2v(:,nyfin-1)
140!  G1%e2f(:,nyfin)=G1%e2f(:,nyfin-1)
141
142
143  CALL agrif_interp(G0%e1t,G1%e1t,'T')   ;   G1%e1t = G1%e1t / REAL(irafx)
144  CALL agrif_interp(G0%e2t,G1%e2t,'T')   ;   G1%e2t = G1%e2t / REAL(irafy)
145
146  CALL agrif_interp(G0%e1u,G1%e1u,'U')   ;   G1%e1u = G1%e1u / REAL(irafx)
147  CALL agrif_interp(G0%e2u,G1%e2u,'U')   ;   G1%e2u = G1%e2u / REAL(irafy) 
148
149  CALL agrif_interp(G0%e1v,G1%e1v,'V')   ;   G1%e1v = G1%e1v / REAL(irafx)
150  CALL agrif_interp(G0%e2v,G1%e2v,'V')   ;   G1%e2v = G1%e2v / REAL(irafy) 
151
152  CALL agrif_interp(G0%e1f,G1%e1f,'F')   ;   G1%e1f = G1%e1f / REAL(irafx)
153  CALL agrif_interp(G0%e2f,G1%e2f,'F')   ;   G1%e2f = G1%e2f / REAL(irafy)
154  !
155  WHERE ( G1%glamt > 180 )   G1%glamt = G1%glamt - 360.
156  WHERE ( G1%glamf > 180 )   G1%glamf = G1%glamf - 360.
157  WHERE ( G1%glamu > 180 )   G1%glamu = G1%glamu - 360.
158  WHERE ( G1%glamv > 180 )   G1%glamv = G1%glamv - 360.
159  !
160  G1%nav_lon=G1%glamt
161  G1%nav_lat=G1%gphit 
162  !             
163  ! Write interpolation result in child coodinates file
164  status = Write_Coordinates(TRIM(Child_filename),G1)
165  !
166  WRITE(*,*) 'Child domain position : '
167  IF( ln_agrif_domain ) THEN
168     WRITE(*,*) 'latmin =',G1%gphit(3,3)
169     WRITE(*,*) 'latmax =',G1%gphit(nxfin-2,nyfin-2)
170     WRITE(*,*) 'lonmin =',G1%glamt(3,3)
171     WRITE(*,*) 'lonmax =',G1%glamt(nxfin-2,nyfin-2)
172  ELSE
173     WRITE(*,*) 'latmin =',G1%gphit(1,1)
174     WRITE(*,*) 'latmax =',G1%gphit(nxfin,nyfin)
175     WRITE(*,*) 'lonmin =',G1%glamt(1,1)
176     WRITE(*,*) 'lonmax =',G1%glamt(nxfin,nyfin)
177  ENDIF
178  STOP
179END PROGRAM create_coordinate
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
Note: See TracBrowser for help on using the repository browser.