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

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

ice restart should work in the nesting tools now. However ocean restart has been broken for some time

File size: 12.7 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  !
83  WRITE(*,*) ''
84  WRITE(*,*) 'Interpolation of restart file : ',TRIM(restart_ice_file)
85  WRITE(*,*) ''
86  !
87  CALL Read_Ncdf_VarName(restart_ice_file,Ncdf_varname)
88  !       
89  CALL set_child_name(parent_coordinate_file,Childcoordinates)   
90  IF( TRIM(parent_bathy_level) /= '' )   CALL set_child_name(parent_bathy_level,Child_Bathy_Level) 
91  IF( TRIM(parent_bathy_meter) /= '' )   CALL set_child_name(parent_bathy_meter,Child_Bathy_Meter)   
92  !
93  ! create this file
94  !
95  CALL set_child_name(restart_ice_file,Child_file)
96  status = nf90_create(Child_file,NF90_WRITE,ncid)
97  status = nf90_close(ncid)
98  WRITE(*,*) 'Child grid restart file name = ',TRIM(Child_file)     
99  WRITE(*,*) ''
100  !
101  ! read dimensions in parent restart file
102  !
103  CALL Read_Ncdf_dim('x',restart_ice_file,x)
104  CALL Read_Ncdf_dim('y',restart_ice_file,y) 
105  CALL Read_Ncdf_dim('numcat',restart_ice_file,z)
106
107  !
108  ! mask initialization for extrapolation and interpolation
109  !
110  WRITE(*,*) 'mask initialisation on coarse and fine grids'
111  !           
112  status = Read_Local_Coordinates(parent_coordinate_file,G0,(/jpizoom,jpjzoom/),(/x,y/))
113  status = Read_Coordinates(Childcoordinates,G1,Pacifique)
114  !
115  !longitude modification if child domain covers Pacific ocean area
116  !     
117  IF( Pacifique ) THEN
118     WHERE( G0%nav_lon < 0 )
119        G0%nav_lon = G0%nav_lon + 360.
120     END WHERE
121     WHERE( G1%nav_lon < 0 )
122        G1%nav_lon = G1%nav_lon + 360.
123     END WHERE
124  ENDIF
125  !
126  ! one needs bathy_level
127  IF( TRIM(parent_bathy_level) /= '' ) THEN
128     status = Read_bathy_level(TRIM(parent_bathy_level),G0)
129     status = Read_bathy_level(TRIM(child_bathy_level),G1)
130  ELSE
131     status = read_bathy_meter(TRIM(parent_bathy_meter),G0)
132     status = read_bathy_meter(TRIM(child_bathy_meter),G1)
133     CALL meter_to_levels(G0)
134     CALL meter_to_levels(G1)
135  ENDIF
136  ! get masks
137  CALL Init_mask(parent_bathy_level,G0,x,y)
138  CALL Init_mask(child_bathy_level,G1,nxfin,nyfin)
139 
140  !
141  ! write dimensions in output file
142  WRITE(*,*) 'write dimensions'
143  !         
144  CALL Write_Ncdf_dim('x',Child_file,nxfin)
145  CALL Write_Ncdf_dim('y',Child_file,nyfin)
146  CALL Write_Ncdf_dim('numcat',Child_file,z)
147  CALL Write_Ncdf_dim('time_counter',Child_file,0) 
148  !
149  !
150  DO ji = 1,SIZE(Ncdf_varname)     
151     !     
152     ! loop on variables names
153     varname = TRIM(Ncdf_varname(ji))
154     !     
155     WRITE(*,*) 'var = ',TRIM(varname)     
156     !
157     SELECT CASE (varname)
158        !
159        !copy nav_lon from child coordinates to output file     
160        !
161     CASE('nav_lon')
162        CALL Read_Ncdf_var(varname,TRIM(Childcoordinates),tabtemp2D) 
163        CALL Write_Ncdf_var(varname,(/'x','y'/),Child_file,tabtemp2D,'float')
164        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file, MINVAL(tabtemp2D),MAXVAL(tabtemp2D))
165        DEALLOCATE(tabtemp2D)
166        Interpolation = .FALSE.
167        !       
168        !copy nav_lat from child coordinates to output file
169        !
170     CASE('nav_lat')             
171        CALL Read_Ncdf_var(varname,TRIM(Childcoordinates),tabtemp2D) 
172        CALL Write_Ncdf_var(varname,(/'x','y'/),Child_file,tabtemp2D,'float')
173        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file,MINVAL(tabtemp2D),MAXVAL(tabtemp2D)) 
174        DEALLOCATE(tabtemp2D)
175        Interpolation = .FALSE.
176        !
177        !copy numcat from restart_ice_file to output file
178        !
179     CASE('numcat')
180        CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),nav_lev) 
181        CALL Write_Ncdf_var(varname,varname,Child_file,nav_lev,'float')
182        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file)     
183        Interpolation = .FALSE.
184        !
185        !copy time from restart_ice_file to output file                       
186        !
187     CASE('time_counter')
188        CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),tabtemp1D) 
189        tabtemp1D = tabtemp1D * rhot
190        CALL Write_Ncdf_var(varname,'time_counter',Child_file,tabtemp1D,'double')
191        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file) 
192        DEALLOCATE(tabtemp1D)
193        Interpolation = .FALSE.
194        !
195        !copy info from restart_ice_file to output file
196        !
197     CASE('kt_ice') 
198        CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),tabtemp0dreal) 
199        tabtemp0dreal = tabtemp0dreal * rhot
200        CALL Write_Ncdf_var(varname,Child_file,tabtemp0dreal,'double')
201        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file) 
202        Interpolation = .FALSE.
203        !
204     CASE('nn_fsbc') 
205        CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),tabtemp0dreal) 
206        CALL Write_Ncdf_var(varname,Child_file,tabtemp0dreal,'double')
207        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file) 
208        Interpolation = .FALSE.
209        !
210     CASE('frc_voltop','frc_volbot','frc_temtop','frc_tembot') 
211        CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),tabtemp0dreal) 
212        CALL Write_Ncdf_var(varname,Child_file,tabtemp0dreal,'double')
213        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file) 
214        Interpolation = .FALSE.
215        !
216     CASE('u_ice')
217        vert_coord_name = '1'             
218        posvar='U'
219        Interpolation = .TRUE.           
220       
221     CASE('v_ice')
222        vert_coord_name = '1'             
223        posvar='V'
224        Interpolation = .TRUE.           
225
226     CASE('stress12_i')
227        vert_coord_name = '1'             
228        posvar='F'
229        Interpolation = .TRUE.           
230       
231     CASE DEFAULT
232        IF( Get_NbDims(TRIM(Ncdf_varname(ji)),TRIM(restart_ice_file)) == 4 ) THEN
233           vert_coord_name = 'numcat'
234        ELSEIF( Get_NbDims(TRIM(Ncdf_varname(ji)),TRIM(restart_ice_file)) == 3 ) THEN
235           vert_coord_name = '1'
236        ENDIF
237        posvar='T'
238        Interpolation = .TRUE.           
239        !     
240     END SELECT
241
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(*,*) '******* restart file successfully created *******' 
342  WRITE(*,*) ' '
343  !
344  STOP
345END PROGRAM create_rstrt_ice
346
347
Note: See TracBrowser for help on using the repository browser.