New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_create_restart_ice.f90 in utils/tools/NESTING/src – NEMO

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

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