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

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

BGC restart should now work in the nesting tool. All of that must of course be tested in realistic simulations. Waiting for the feedbacks…

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