source: utils/tools/NESTING/src/agrif_create_restart_ice.f90 @ 10393

Last change on this file since 10393 was 10393, checked in by clem, 23 months ago

ocean restart from nesting tool should now work. BGC restart must still be debugged but it is out of my expertise

File size: 12.9 KB
Line 
1!
2!************************************************************************
3! Fortran 95 OPA Nesting tools                  *
4!                          *
5!     Copyright (C) 2005 Florian Lemari�(Florian.Lemarie@imag.fr) *
6!                          *
7!************************************************************************
8!
9PROGRAM create_rstrt_ice
10  !
11  USE NETCDF
12  USE bilinear_interp
13  USE bicubic_interp
14  USE agrif_readwrite
15  USE io_netcdf 
16  USE agrif_extrapolation
17  USE agrif_interpolation
18  USE agrif_partial_steps
19  USE agrif_connect_topo
20  !
21  IMPLICIT NONE
22  !
23  !************************************************************************
24  !                           *
25  ! PROGRAM  CREATE_RSTRT_ICE                   *
26  !                           *
27  ! program to interpolate parent grid restart file to child grid    *
28  !                           *
29  !                           *
30  !Interpolation is carried out using bilinear interpolation      *
31  !routine from SCRIP package                *     
32  !                           *
33  !http://climate.lanl.gov/Software/SCRIP/            *
34  !************************************************************************
35  !
36  ! variables declaration
37  !     
38  CHARACTER*20,DIMENSION(:),POINTER :: Ncdf_varname => NULL()
39  CHARACTER*20 :: vert_coord_name
40  CHARACTER*1 :: posvar
41  CHARACTER*100 :: Child_file,Childcoordinates,varname,Child_Bathy_Level,Child_Bathy_Meter   
42  REAL*8, POINTER, DIMENSION(:,:,:) :: tabvar00, tabvar3d,tabinterp3d,mask => NULL()
43  REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabinterp4d,tabvar0,tabvar1,tabvar2,tabvar3 => NULL()
44  REAL*8, POINTER, DIMENSION(:) :: tabtemp1D,nav_lev => NULL()
45  INTEGER, POINTER, DIMENSION(:) :: tabtemp1DInt => NULL()
46  REAL*8, POINTER, DIMENSION(:,:) :: tabtemp2D,zwf => NULL()
47  REAL*8, POINTER, DIMENSION(:,:) :: e1f,e2f,e1v,e2v,e1u,e2u,e1t,e2t => NULL()
48  REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabtemp4D => NULL()
49  INTEGER,DIMENSION(:),POINTER :: src_add,dst_add  => NULL()
50  REAL*8,DIMENSION(:,:),POINTER :: matrix => NULL()
51  LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL()
52  LOGICAL, DIMENSION(:,:,:), POINTER :: detected_pts
53  LOGICAL :: Interpolation,Extrapolation,Pacifique
54  INTEGER :: narg,iargc,ncid,x,y,z
55  INTEGER :: ji,jl,status,varid,numdims
56  CHARACTER(len=20),DIMENSION(4) :: dimnames
57  CHARACTER(len=80) :: namelistname
58  TYPE(Coordinates) :: G0,G1
59  INTEGER :: jpl,jpiglo,jpjglo
60  REAL*8 :: tabtemp0dreal
61  INTEGER :: tabtemp0dint
62
63  LOGICAL, PARAMETER :: conservation = .FALSE.
64  !       
65  narg = iargc()
66  IF (narg == 0) THEN
67     namelistname = 'namelist.input'
68  ELSE
69     CALL getarg(1,namelistname)
70  ENDIF
71  !
72  ! read input file
73  !
74  CALL read_namelist(namelistname)
75  !     
76  IF(TRIM(restart_ice_file) == '') THEN
77     WRITE(*,*) 'no ice restart file specified in ',TRIM(namelistname)
78     STOP
79  ENDIF
80
81  !
82  WRITE(*,*) ''
83  WRITE(*,*) 'Interpolation of restart file : ',TRIM(restart_ice_file)
84  WRITE(*,*) ''
85  !
86  CALL Read_Ncdf_VarName(restart_ice_file,Ncdf_varname)
87  !       
88  CALL set_child_name(parent_coordinate_file,Childcoordinates)   
89  IF( TRIM(parent_bathy_level) /= '' )   CALL set_child_name(parent_bathy_level,Child_Bathy_Level) 
90  IF( TRIM(parent_bathy_meter) /= '' )   CALL set_child_name(parent_bathy_meter,Child_Bathy_Meter)   
91  !
92  ! create this file
93  !
94  CALL set_child_name(restart_ice_file,Child_file)
95  status = nf90_create(Child_file,NF90_WRITE,ncid)
96  status = nf90_close(ncid)
97  WRITE(*,*) 'Child grid restart file name = ',TRIM(Child_file)     
98  WRITE(*,*) ''
99  !
100  ! read dimensions in parent restart file
101  !
102  CALL Read_Ncdf_dim('x',restart_ice_file,x)
103  CALL Read_Ncdf_dim('y',restart_ice_file,y) 
104  CALL Read_Ncdf_dim('numcat',restart_ice_file,z)
105
106  !
107  ! mask initialization for extrapolation and interpolation
108  !
109  WRITE(*,*) 'mask initialisation on coarse and fine grids'
110  !           
111  status = Read_Local_Coordinates(parent_coordinate_file,G0,(/jpizoom,jpjzoom/),(/x,y/))
112  status = Read_Coordinates(Childcoordinates,G1,Pacifique)
113  !
114  !longitude modification if child domain covers Pacific ocean area
115  !     
116  IF( Pacifique ) THEN
117     WHERE( G0%nav_lon < 0 )
118        G0%nav_lon = G0%nav_lon + 360.
119     END WHERE
120     WHERE( G1%nav_lon < 0 )
121        G1%nav_lon = G1%nav_lon + 360.
122     END WHERE
123  ENDIF
124  !
125  ! one needs bathy_level
126  IF( TRIM(parent_bathy_level) /= '' ) THEN
127     status = Read_bathy_level(TRIM(parent_bathy_level),G0)
128     status = Read_bathy_level(TRIM(child_bathy_level),G1)
129  ELSE
130     status = read_bathy_meter(TRIM(parent_bathy_meter),G0)
131     status = read_bathy_meter(TRIM(child_bathy_meter),G1)
132     CALL meter_to_levels(G0)
133     CALL meter_to_levels(G1)
134  ENDIF
135  ! get masks
136  CALL Init_mask(parent_bathy_level,G0,x,y)
137  CALL Init_mask(child_bathy_level,G1,nxfin,nyfin)
138 
139  !
140  ! write dimensions in output file
141  WRITE(*,*) 'write dimensions'
142  !         
143  CALL Write_Ncdf_dim('x',Child_file,nxfin)
144  CALL Write_Ncdf_dim('y',Child_file,nyfin)
145  CALL Write_Ncdf_dim('numcat',Child_file,z)
146  CALL Write_Ncdf_dim('time_counter',Child_file,0) 
147  !
148  !
149  DO ji = 1,SIZE(Ncdf_varname)     
150     !     
151     ! loop on variables names
152     varname = TRIM(Ncdf_varname(ji))
153     !     
154     WRITE(*,*) 'var = ',TRIM(varname)     
155     !
156     SELECT CASE (varname)
157        !
158        !copy nav_lon from child coordinates to output file     
159        !
160     CASE('nav_lon')
161        CALL Read_Ncdf_var(varname,TRIM(Childcoordinates),tabtemp2D) 
162        CALL Write_Ncdf_var(varname,(/'x','y'/),Child_file,tabtemp2D,'float')
163        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file, MINVAL(tabtemp2D),MAXVAL(tabtemp2D))
164        DEALLOCATE(tabtemp2D)
165        Interpolation = .FALSE.
166        !       
167        !copy nav_lat from child coordinates to output file
168        !
169     CASE('nav_lat')             
170        CALL Read_Ncdf_var(varname,TRIM(Childcoordinates),tabtemp2D) 
171        CALL Write_Ncdf_var(varname,(/'x','y'/),Child_file,tabtemp2D,'float')
172        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file,MINVAL(tabtemp2D),MAXVAL(tabtemp2D)) 
173        DEALLOCATE(tabtemp2D)
174        Interpolation = .FALSE.
175        !
176        !copy numcat from restart_ice_file to output file
177        !
178     CASE('numcat')
179        CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),nav_lev) 
180        CALL Write_Ncdf_var(varname,varname,Child_file,nav_lev,'float')
181        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file)     
182        Interpolation = .FALSE.
183        !
184        !copy time from restart_ice_file to output file                       
185        !
186     CASE('time_counter')
187        CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),tabtemp1D) 
188        tabtemp1D = tabtemp1D * rhot
189        CALL Write_Ncdf_var(varname,'time_counter',Child_file,tabtemp1D,'double')
190        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file) 
191        DEALLOCATE(tabtemp1D)
192        Interpolation = .FALSE.
193        !
194        !copy info from restart_ice_file to output file
195        !
196     CASE('kt_ice') 
197        CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),tabtemp0dreal) 
198        tabtemp0dreal = tabtemp0dreal * rhot
199        CALL Write_Ncdf_var(varname,Child_file,tabtemp0dreal,'double')
200        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file) 
201        Interpolation = .FALSE.
202        !
203     CASE('nn_fsbc') 
204        CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),tabtemp0dreal) 
205        CALL Write_Ncdf_var(varname,Child_file,tabtemp0dreal,'double')
206        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file) 
207        Interpolation = .FALSE.
208        !
209     CASE('frc_voltop','frc_volbot','frc_temtop','frc_tembot') 
210        CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),tabtemp0dreal) 
211        CALL Write_Ncdf_var(varname,Child_file,tabtemp0dreal,'double')
212        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file) 
213        Interpolation = .FALSE.
214        !
215     CASE('u_ice')
216        vert_coord_name = '1'             
217        posvar='U'
218        Interpolation = .TRUE.           
219       
220     CASE('v_ice')
221        vert_coord_name = '1'             
222        posvar='V'
223        Interpolation = .TRUE.           
224
225     CASE('stress12_i')
226        vert_coord_name = '1'             
227        posvar='F'
228        Interpolation = .TRUE.           
229       
230     CASE DEFAULT
231        IF( Get_NbDims(TRIM(varname),TRIM(restart_ice_file)) == 4 ) THEN
232           vert_coord_name = 'numcat'
233        ELSEIF( Get_NbDims(TRIM(varname),TRIM(restart_ice_file)) == 3 ) THEN
234           vert_coord_name = '1'
235        ENDIF
236        posvar='T'
237        Interpolation = .TRUE.           
238        !     
239     END SELECT
240
241     ! --- start interpolation --- !
242     IF( Interpolation ) THEN
243        !
244        !
245        IF( vert_coord_name == '1' ) THEN
246           jpl = 1
247        ELSE
248           jpl = z
249        ENDIF
250        !
251        ALLOCATE(detected_pts(SIZE(G0%tmask,1),SIZE(G0%tmask,2),jpl))                             
252        ALLOCATE(tabvar00(x,y,1))
253        ALLOCATE(tabvar0(x,y,jpl,1))
254        ALLOCATE(tabvar1(x,y,1,1))
255        ALLOCATE(tabvar2(x,y,1,1))
256        ALLOCATE(tabvar3(x,y,1,1))
257        ALLOCATE(masksrc(x,y))
258        ALLOCATE(tabinterp4d(nxfin,nyfin,jpl,1)) 
259
260        IF( vert_coord_name == '1' ) THEN
261           CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),tabvar00)
262        ELSE
263           CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),tabvar0)
264        ENDIF
265
266        DO jl = 1, jpl
267           !
268           WRITE(*,*) 'interpolate/extrapolate for category = ',jl   
269           !
270           IF( vert_coord_name == '1' ) THEN
271              tabvar1(:,:,1,1) = tabvar00(:,:,1)
272              tabvar2(:,:,1,1) = tabvar00(:,:,1)
273              tabvar3(:,:,1,1) = tabvar00(:,:,1)
274           ELSE
275              tabvar1(:,:,1,1) = tabvar0(:,:,jl,1)
276              tabvar2(:,:,1,1) = tabvar0(:,:,jl,1)
277              tabvar3(:,:,1,1) = tabvar0(:,:,jl,1)
278           ENDIF
279           !
280           CALL extrap_detect(G0,G1,detected_pts(:,:,jl),1,posvar)                                           
281           CALL correct_field(detected_pts(:,:,jl),tabvar1,tabvar2,tabvar3,G0,nav_lev,masksrc,1,posvar) !clem: nav_lev is useless here
282           
283           SELECT CASE(TRIM(interp_type))
284           CASE('bilinear')                                                       
285              CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
286              CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,jl,1),nxfin,nyfin,matrix,src_add,dst_add)
287           CASE('bicubic')                                   
288              CALL get_remap_bicub(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
289              CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,jl,1),nxfin,nyfin,matrix,src_add,dst_add) 
290           END SELECT
291           !
292           IF( conservation ) THEN ! clem: currently this does not work (and is only coded for T)
293              SELECT CASE (posvar)
294              CASE( 'T' )
295                 CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,jl,1), &
296                    &                        G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
297              CASE( 'U' )
298                 CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,jl,1), &
299                    &                        G0%e1u,G0%e2u,G1%e1u,G1%e2u,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
300              CASE( 'V' )
301                 CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,jl,1), &
302                    &                        G0%e1v,G0%e2v,G1%e1v,G1%e2v,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
303              CASE( 'F' )
304                 CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,jl,1), &
305                    &                        G0%e1f,G0%e2f,G1%e1f,G1%e2f,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
306              END SELECT
307           ENDIF
308           tabinterp4d(:,:,jl,1) =  tabinterp4d(:,:,jl,1) * G1%tmask(:,:,1)
309
310           IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,src_add,dst_add)   
311
312        ENDDO
313
314        dimnames(1)='x'
315        dimnames(2)='y'
316        IF( vert_coord_name == '1' ) THEN
317           dimnames(3)='time_counter'
318
319           ALLOCATE(tabvar3d(SIZE(tabinterp4d,1),SIZE(tabinterp4d,2),SIZE(tabinterp4d,3)))
320           tabvar3d=tabinterp4d(:,:,:,1)
321           CALL Write_Ncdf_var(varname,dimnames,Child_file,tabvar3d,'double')
322           DEALLOCATE(tabvar3d)
323        ELSE
324           dimnames(3)='numcat'
325           dimnames(4)='time_counter'
326           
327           CALL Write_Ncdf_var(TRIM(varname),dimnames,Child_file,tabinterp4d,'double')
328        ENDIF
329        !             
330        CALL Copy_Ncdf_att(TRIM(varname),TRIM(restart_ice_file),Child_file)
331
332        DEALLOCATE(detected_pts)
333        DEALLOCATE(tabinterp4d)
334        DEALLOCATE(tabvar00, tabvar0, tabvar1,tabvar2,tabvar3)
335        DEALLOCATE(masksrc) 
336       
337     ENDIF
338  END DO
339  !
340  WRITE(*,*) ' '
341  WRITE(*,*) ' --- list of all the variables that have been interpolated --- '
342  WRITE(*,*) Ncdf_varname
343  WRITE(*,*) ' '
344  WRITE(*,*) '******* restart file successfully created *******' 
345  WRITE(*,*) ' '
346  !
347  STOP
348END PROGRAM create_rstrt_ice
349
350
Note: See TracBrowser for help on using the repository browser.