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

Last change on this file since 12253 was 10398, checked in by clem, 2 years 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.