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 @ 10393

Last change on this file since 10393 was 10393, checked in by clem, 5 years 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
RevLine 
[10383]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
[10393]231        IF( Get_NbDims(TRIM(varname),TRIM(restart_ice_file)) == 4 ) THEN
[10383]232           vert_coord_name = 'numcat'
[10393]233        ELSEIF( Get_NbDims(TRIM(varname),TRIM(restart_ice_file)) == 3 ) THEN
[10383]234           vert_coord_name = '1'
235        ENDIF
236        posvar='T'
237        Interpolation = .TRUE.           
238        !     
239     END SELECT
240
[10393]241     ! --- start interpolation --- !
[10383]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(*,*) ' '
[10393]341  WRITE(*,*) ' --- list of all the variables that have been interpolated --- '
342  WRITE(*,*) Ncdf_varname
343  WRITE(*,*) ' '
[10383]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.