source: utils/tools/NESTING/src/agrif_create_restart.f90 @ 10393

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

ocean restart from nesting tool should now work. BGC restart must still be debugged but it is out of my expertise

  • Property svn:keywords set to Id
File size: 21.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
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                    *
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(:,:) :: lonChild,latChild => NULL()
43  REAL*8, POINTER, DIMENSION(:,:) :: lonParent,latParent => NULL()
44  REAL*8, POINTER, DIMENSION(:,:,:) :: tabvar3d,tabinterp3d,mask => NULL()
45  REAL*8, POINTER, DIMENSION(:,:,:) :: fmask,fse3u,fse3v,fse3t => NULL()
46  REAL*8, POINTER, DIMENSION(:,:,:,:) :: un,ub,vn,vb,tn,tb,sn,sb,e3t_n,e3t_b => NULL()
47  REAL*8, POINTER, DIMENSION(:,:,:) :: sshn,sshb => NULL()     
48  REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabinterp4d,tabvar1,tabvar2,tabvar3 => NULL()
49  REAL*8, POINTER, DIMENSION(:) :: timedepth_temp => NULL()
50  REAL*8, POINTER, DIMENSION(:) :: tabtemp1D,nav_lev => NULL()
51  INTEGER, POINTER, DIMENSION(:) :: tabtemp1DInt => NULL()
52  REAL*8, POINTER, DIMENSION(:,:) :: tabtemp2D,zwf => NULL()
53  REAL*8, POINTER, DIMENSION(:,:) :: e1f,e2f,e1v,e2v,e1u,e2u,e1t,e2t => NULL()
54  REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabtemp4D => NULL()
55  INTEGER,DIMENSION(:),POINTER :: src_add,dst_add  => NULL()
56  REAL*8,DIMENSION(:,:),POINTER :: matrix,bathy_G0 => NULL()
57  LOGICAL,DIMENSION(:,:,:),POINTER :: Tmask => NULL()
58  LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL()
59  LOGICAL, DIMENSION(:,:,:), POINTER :: detected_pts
60  LOGICAL :: Interpolation,Extrapolation,Pacifique,op
61  REAL*8 :: za1,za0,zsur,zacr,zkth,zdepth,zdepwp,zmin,zmax,zdiff,ze3tp,ze3wp       
62  INTEGER :: narg,iargc,ncid,x,y,z,z_a,x_a,y_a,z_b,nbvert_lev,m
63  REAL*8 :: now_wght,before_wght
64  INTEGER :: i,j,k,status,ji,jj
65  CHARACTER(len=20),DIMENSION(4) :: dimnames
66  CHARACTER(len=80) :: namelistname
67  TYPE(Coordinates) :: G0,G1
68  INTEGER :: jpi,jpj,jpk,jpni,jpnj,jpnij,jpiglo,jpjglo,nlcit,nlcjt,nldit
69  INTEGER :: nldjt,nleit,nlejt,nimppt,njmppt
70  REAL*8 :: tabtemp0dreal
71  INTEGER :: tabtemp0dint
72  CHARACTER(len=20) :: timedimname
73
74  LOGICAL, PARAMETER :: conservation = .FALSE.
75  !     
76  ! Variables for dimg
77  !
78!  INTEGER :: ino0, it0, ipcg0, isor0, itke0,nfice,nfbulk
79!  INTEGER :: irecl8, irec,ndastp,narea,nsolv
80!  INTEGER :: jk,kt            ! dummy loop indices
81!  INTEGER :: inum             ! temporary logical unit
82!  INTEGER :: ios1 , ios2      ! flag for ice and bulk in the current run
83!  INTEGER :: ios3             ! flag for free surface.  0 = none 1 = yes.  0 = none 1 = yes
84!  INTEGER :: ios4             ! flag for coupled (1) or not (0)   
85  !       
86  narg = iargc()
87  IF (narg == 0) THEN
88     namelistname = 'namelist.input'
89  ELSE
90     CALL getarg(1,namelistname)
91  ENDIF
92  !
93  ! read input file
94  !
95  CALL read_namelist(namelistname)
96  !     
97  IF(TRIM(restart_file) == '') THEN
98     WRITE(*,*) 'no restart file specified in ',TRIM(namelistname)
99     STOP
100  END IF
101
102  IF (iom_activated) THEN
103     timedimname = 'time_counter'
104  ELSE
105     timedimname='time'
106  ENDIF
107
108  !
109  WRITE(*,*) ''
110  WRITE(*,*) 'Interpolation of restart file : ',TRIM(restart_file)
111  WRITE(*,*) ''
112  !
113  CALL Read_Ncdf_VarName(restart_file,Ncdf_varname)
114  !       
115  CALL set_child_name(parent_coordinate_file,Childcoordinates)   
116  IF( TRIM(parent_bathy_level) /= '' )   CALL set_child_name(parent_bathy_level,Child_Bathy_Level) 
117  IF( TRIM(parent_bathy_meter) /= '' )   CALL set_child_name(parent_bathy_meter,Child_Bathy_Meter)   
118  !
119  ! create this file
120  !
121  CALL set_child_name(restart_file,Child_file)
122  status = nf90_create(Child_file,NF90_WRITE,ncid)
123  status = nf90_close(ncid)
124  WRITE(*,*) 'Child grid restart file name = ',TRIM(Child_file)     
125  WRITE(*,*) ''
126
127  !
128  ! read dimensions in parent restart file
129  !
130  CALL Read_Ncdf_dim('x',restart_file,x)
131  CALL Read_Ncdf_dim('y',restart_file,y) 
132  CALL Read_Ncdf_dim('nav_lev',restart_file,z)
133  IF (.NOT.iom_activated) THEN
134     CALL Read_Ncdf_dim('x_a',restart_file,x_a)
135     CALL Read_Ncdf_dim('y_a',restart_file,y_a)
136     CALL Read_Ncdf_dim('z_a',restart_file,z_a)
137     CALL Read_Ncdf_dim('z_b',restart_file,z_b)
138  ENDIF
139
140  IF( z .NE. N ) THEN
141     WRITE(*,*) '***'
142     WRITE(*,*) 'Number of vertical levels doesn t match between namelist and restart file'
143     WRITE(*,*) 'Please check the values in namelist file'
144     STOP
145  ENDIF
146  !
147  ! mask initialization for extrapolation and interpolation
148  !
149  WRITE(*,*) 'mask initialisation on coarse and fine grids'
150  !           
151  status = Read_Local_Coordinates(parent_coordinate_file,G0,(/jpizoom,jpjzoom/),(/x,y/))
152  status = Read_Coordinates(Childcoordinates,G1,Pacifique)
153  !
154  !longitude modification if child domain covers Pacific ocean area
155  !     
156  IF( Pacifique ) THEN
157     !
158     WHERE( G0%nav_lon < 0 )
159        G0%nav_lon = G0%nav_lon + 360.
160     END WHERE
161     !             
162     WHERE( G1%nav_lon < 0 )
163        G1%nav_lon = G1%nav_lon + 360.
164     END WHERE
165     !
166  ENDIF
167  !
168  ! one needs bathy_level
169  IF( TRIM(parent_bathy_level) /= '' ) THEN
170     status = Read_bathy_level(TRIM(parent_bathy_level),G0)
171     status = Read_bathy_level(TRIM(child_bathy_level),G1)
172  ELSE
173     status = read_bathy_meter(TRIM(parent_bathy_meter),G0)
174     status = read_bathy_meter(TRIM(child_bathy_meter),G1)
175     CALL meter_to_levels(G0)
176     CALL meter_to_levels(G1)
177  ENDIF
178  ! get masks
179  CALL Init_mask(parent_bathy_level,G0,x,y)
180  CALL Init_mask(child_bathy_level,G1,1,1)
181
182  G0%tmask = 1.   
183
184  DO k=1,z
185     ALLOCATE(tabvar1(x,y,1,1))
186     CALL Read_Ncdf_var('sn',TRIM(restart_file),tabvar1,1,k)
187     WHERE( tabvar1(:,:,1,1) == 0. ) 
188        G0%tmask(:,:,k) = 0.
189     END WHERE
190     DEALLOCATE(tabvar1)
191  END DO
192  !
193  G0%umask(1:x-1,:,:) = G0%tmask(1:x-1,:,:)*G0%tmask(2:x,:,:)
194  G0%umask(x,:,:)     = G0%tmask(x,:,:)
195  G0%vmask(:,1:y-1,:) = G0%tmask(:,1:y-1,:)*G0%tmask(:,2:y,:)
196  G0%vmask(:,y,:)     = G0%tmask(:,y,:)
197  !     
198  G0%fmask(1:x-1,1:y-1,:) = G0%tmask(1:x-1,1:y-1,:)*G0%tmask(2:x,1:y-1,:)* &
199     &                      G0%tmask(1:x-1,2:y,:)*G0%tmask(2:x,2:y,:) 
200  G0%fmask(x,:,:) = G0%tmask(x,:,:)
201  G0%fmask(:,y,:) = G0%tmask(:,y,:)
202  !
203  !
204  ! write dimensions in output file
205  WRITE(*,*) 'write dimensions'
206  !         
207  CALL Write_Ncdf_dim('x',Child_file,nxfin)
208  CALL Write_Ncdf_dim('y',Child_file,nyfin)
209  CALL Write_Ncdf_dim('nav_lev',Child_file,z)
210  CALL Write_Ncdf_dim(TRIM(timedimname),Child_file,0) 
211  IF (.NOT.iom_activated) THEN
212     CALL Write_Ncdf_dim('x_a',Child_file,x_a)
213     CALL Write_Ncdf_dim('y_a',Child_file,y_a)
214     CALL Write_Ncdf_dim('z_a',Child_file,z_a)
215     CALL Write_Ncdf_dim('z_b',Child_file,z_b)
216  ENDIF
217  !
218  !
219  !
220  !
221  DO ji = 1,SIZE(Ncdf_varname)     
222     !     
223     ! loop on variables names
224     varname = TRIM(Ncdf_varname(ji))
225     WRITE(*,*) 'var = ',TRIM(varname)     
226     !     
227     SELECT CASE (TRIM(varname))
228        !
229     CASE('nav_lon')
230        CALL Read_Ncdf_var('nav_lon',TRIM(Childcoordinates),tabtemp2D) 
231        CALL Write_Ncdf_var('nav_lon',(/'x','y'/),Child_file,tabtemp2D,'float')
232        CALL Copy_Ncdf_att('nav_lon',TRIM(restart_file),Child_file,MINVAL(tabtemp2D),MAXVAL(tabtemp2D))
233        DEALLOCATE(tabtemp2D)
234        Interpolation = .FALSE.
235        !       
236     CASE('nav_lat')             
237        CALL Read_Ncdf_var('nav_lat',TRIM(Childcoordinates),tabtemp2D) 
238        CALL Write_Ncdf_var('nav_lat',(/'x','y'/),Child_file,tabtemp2D,'float')
239        CALL Copy_Ncdf_att('nav_lat',TRIM(restart_file),Child_file,MINVAL(tabtemp2D),MAXVAL(tabtemp2D)) 
240        DEALLOCATE(tabtemp2D)
241        Interpolation = .FALSE.
242        !
243     CASE('nav_lev')
244        CALL Read_Ncdf_var('nav_lev',TRIM(restart_file),nav_lev) 
245        CALL Write_Ncdf_var('nav_lev','nav_lev',Child_file,nav_lev,'float')
246        CALL Copy_Ncdf_att('nav_lev',TRIM(restart_file),Child_file)     
247        Interpolation = .FALSE.
248        !
249     CASE('time')
250        CALL Read_Ncdf_var('time',TRIM(restart_file),tabtemp1D) 
251        CALL Write_Ncdf_var('time',TRIM(timedimname),Child_file,tabtemp1D,'float')
252        CALL Copy_Ncdf_att('time',TRIM(restart_file),Child_file) 
253        DEALLOCATE(tabtemp1D)
254        Interpolation = .FALSE.
255        !
256     CASE('time_counter')
257        CALL Read_Ncdf_var('time_counter',TRIM(restart_file),tabtemp1D) 
258        tabtemp1D = tabtemp1D * rhot
259        CALL Write_Ncdf_var('time_counter',TRIM(timedimname),Child_file,tabtemp1D,'double')
260        CALL Copy_Ncdf_att('time_counter',TRIM(restart_file),Child_file) 
261        DEALLOCATE(tabtemp1D)
262        Interpolation = .FALSE.
263        !
264     CASE('kt','ndastp','adatrj','ntime','nn_fsbc','rdt') 
265        IF (iom_activated) THEN
266           CALL Read_Ncdf_var(TRIM(varname),TRIM(restart_file),tabtemp0dreal) 
267           SELECT CASE (TRIM(varname))
268           CASE('rdt')
269              tabtemp0dreal = tabtemp0dreal / rhot
270           CASE('kt')
271              tabtemp0dreal = tabtemp0dreal * rhot
272           END SELECT
273           CALL Write_Ncdf_var(TRIM(varname),Child_file,tabtemp0dreal,'double')
274        ELSE
275           CALL Read_Ncdf_var(TRIM(varname),TRIM(restart_file),tabtemp4D) 
276           dimnames(1)='x_a'
277           dimnames(2)='y_a'
278           dimnames(3)='z_b'
279           dimnames(4)=TRIM(timedimname)           
280           CALL Write_Ncdf_var(TRIM(varname),dimnames,Child_file,tabtemp4D,'double')
281           DEALLOCATE(tabtemp4D)
282        ENDIF
283        CALL Copy_Ncdf_att(TRIM(varname),TRIM(restart_file),Child_file) 
284        Interpolation = .FALSE.
285        !
286     CASE('frc_v','frc_t','frc_s') 
287        CALL Read_Ncdf_var(varname,TRIM(restart_ice_file),tabtemp0dreal) 
288        CALL Write_Ncdf_var(varname,Child_file,tabtemp0dreal,'double')
289        CALL Copy_Ncdf_att(varname,TRIM(restart_ice_file),Child_file) 
290        Interpolation = .FALSE.
291        !
292        ! Variable interpolation according to their position on grid
293        !                                   
294     CASE('ssu_m','utau_b','un_bf','un','ub') 
295        IF( Get_NbDims(TRIM(varname),TRIM(restart_file)) == 4 ) THEN
296           vert_coord_name = 'nav_lev'
297        ELSEIF( Get_NbDims(TRIM(varname),TRIM(restart_file)) == 3 ) THEN
298           vert_coord_name = '1'
299        ENDIF
300        posvar='U'
301        Interpolation = .TRUE. 
302        !
303     CASE('ssv_m','vtau_b','vn_bf','vn','vb') 
304        IF( Get_NbDims(TRIM(varname),TRIM(restart_file)) == 4 ) THEN
305           vert_coord_name = 'nav_lev'
306        ELSEIF( Get_NbDims(TRIM(varname),TRIM(restart_file)) == 3 ) THEN
307           vert_coord_name = '1'
308        ENDIF
309        posvar='V'
310        Interpolation = .TRUE. 
311        !
312     CASE DEFAULT
313        IF( Get_NbDims(TRIM(varname),TRIM(restart_file)) == 4 ) THEN
314           vert_coord_name = 'nav_lev'
315        ELSEIF( Get_NbDims(TRIM(varname),TRIM(restart_file)) == 3 ) THEN
316           vert_coord_name = '1'
317        ENDIF
318        posvar='T'
319        Interpolation = .TRUE.           
320        !     
321     END SELECT
322     !     
323     ! --- start interpolation --- !
324     IF( Interpolation ) THEN
325        !       
326        IF( vert_coord_name == '1' ) THEN
327           nbvert_lev = 1
328        ELSE
329           nbvert_lev = z
330        ENDIF
331        !
332
333        ALLOCATE(detected_pts(SIZE(G0%tmask,1),SIZE(G0%tmask,2),nbvert_lev))                             
334        ALLOCATE(tabvar1(x,y,1,2))
335        ALLOCATE(tabvar2(x,y,1,1))
336        ALLOCATE(tabvar3(x,y,1,1))
337        ALLOCATE(masksrc(x,y))
338        ALLOCATE(tabinterp4d(nxfin,nyfin,1,1)) 
339
340        !       
341        DO n = 1,nbvert_lev
342           !
343           WRITE(*,*) 'interpolate/extrapolate for vertical level = ',n   
344           !                           
345           CALL Read_Ncdf_var(varname,TRIM(restart_file),tabvar1,1,n)
346           IF(n==1) THEN
347              !                           
348           ELSE IF (n==2) THEN
349              tabvar2(:,:,:,1) = tabvar1(:,:,:,2)
350           ELSE
351              tabvar3(:,:,:,1) = tabvar2(:,:,:,1)
352              tabvar2(:,:,:,1) = tabvar1(:,:,:,2)
353
354           ENDIF
355           !                           
356           SELECT CASE(posvar)
357              !
358           CASE('T')
359              !
360              IF(MAXVAL(G1%tmask(:,:,n)) == 0.) THEN           
361                 tabinterp4d = 0.0
362                 WRITE(*,*) 'only land points on level ',n                     
363              ELSE                       
364
365                 CALL extrap_detect(G0,G1,detected_pts(:,:,n),n)                                           
366
367                 CALL correct_field(detected_pts(:,:,n),tabvar1,tabvar2,tabvar3,G0,nav_lev,masksrc,n)                                 
368
369                 ! for the following variables, you do not want to mask the values
370                 IF(  TRIM(varname) == 'e3t_n' .OR. TRIM(varname) == 'e3t_b' .OR. TRIM(varname) == 'e3t_m' .OR. &
371                    & TRIM(varname) == 'fraqsr_1lev' .OR. TRIM(varname) == 'frq_m' .OR. TRIM(varname) == 'surf_ini' ) THEN
372                    masksrc(:,:) = .TRUE.
373                 ENDIF
374                 
375                 SELECT CASE(TRIM(interp_type))
376                 CASE('bilinear')                                                       
377                    CALL get_remap_matrix(G0%nav_lat,G1%nav_lat, &
378                       G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
379                    CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
380                       matrix,src_add,dst_add)     
381                 CASE('bicubic')                                   
382                    CALL get_remap_bicub(G0%nav_lat,G1%nav_lat, &
383                       G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
384                    CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),&
385                       nxfin,nyfin,matrix,src_add,dst_add) 
386                 END SELECT
387                 !
388                 IF( conservation ) THEN ! clem: it currently does not work
389                    CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
390                       G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
391                 ENDIF
392
393              ENDIF
394             
395              IF( ALL(masksrc) ) THEN
396                 tabinterp4d(:,:,1,1) =  tabinterp4d(:,:,1,1)
397              ELSE
398                 tabinterp4d(:,:,1,1) =  tabinterp4d(:,:,1,1) * G1%tmask(:,:,n)
399              ENDIF
400              !
401           CASE('U')
402              !
403              IF(MAXVAL(G1%umask(:,:,n)) == 0) THEN
404                 tabinterp4d = 0.0
405                 WRITE(*,*) 'only land points on level ',n 
406              ELSE     
407                 !
408                 CALL extrap_detect(G0,G1,detected_pts(:,:,n),n,'U') 
409                 CALL correct_field(detected_pts(:,:,n),tabvar1,tabvar2,tabvar3,G0,nav_lev,masksrc,n,'U')
410                 !                                                           
411                 SELECT CASE(TRIM(interp_type))
412                 CASE('bilinear')                                                       
413                    CALL get_remap_matrix(G0%gphiu,G1%gphiu,   &
414                       G0%glamu,G1%glamu,masksrc,matrix,src_add,dst_add)
415                    CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
416                       matrix,src_add,dst_add)     
417                 CASE('bicubic')                                   
418                    CALL get_remap_bicub(G0%gphiu,G1%gphiu,   &
419                       G0%glamu,G1%glamu,masksrc,matrix,src_add,dst_add)
420                    CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),&
421                       nxfin,nyfin,matrix,src_add,dst_add)                       
422                 END SELECT
423                 !                     
424                 IF( conservation ) THEN ! clem: not coded for U
425                    CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
426                       G0%e1u,G0%e2u,G1%e1u,G1%e2u,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
427                 ENDIF
428              ENDIF
429
430              tabinterp4d(:,:,1,1) =  tabinterp4d(:,:,1,1) * G1%umask(:,:,n)
431              !
432           CASE('V')
433              !     
434              IF(MAXVAL(G1%vmask(:,:,n)) == 0) THEN
435                 tabinterp4d = 0.0
436                 WRITE(*,*) 'only land points on level ',n 
437              ELSE     
438                 !
439
440                 CALL extrap_detect(G0,G1,detected_pts(:,:,n),n,'V')
441
442                 CALL correct_field(detected_pts(:,:,n),tabvar1,tabvar2,tabvar3,G0,nav_lev,masksrc,n,'V')
443                 !                                                           
444                 SELECT CASE(TRIM(interp_type))
445                 CASE('bilinear')                                                       
446                    CALL get_remap_matrix(G0%gphiv,G1%gphiv,   &
447                       G0%glamv,G1%glamv,masksrc,matrix,src_add,dst_add)
448                    CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
449                       matrix,src_add,dst_add)     
450                 CASE('bicubic')                                   
451                    CALL get_remap_bicub(G0%gphiv,G1%gphiv,   &
452                       G0%glamv,G1%glamv,masksrc,matrix,src_add,dst_add)
453                    CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),&
454                       nxfin,nyfin,matrix,src_add,dst_add)                       
455                 END SELECT
456                 !                     
457                 IF( conservation ) THEN ! clem: not coded for V
458                    CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
459                       G0%e1v,G0%e2v,G1%e1v,G1%e2v,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
460                 ENDIF
461              ENDIF
462
463              tabinterp4d(:,:,1,1) =  tabinterp4d(:,:,1,1) * G1%vmask(:,:,n)
464              !
465           END SELECT
466           !
467           !
468           dimnames(1)='x'
469           dimnames(2)='y'
470           IF( vert_coord_name == '1' ) THEN
471              dimnames(3)=TRIM(timedimname)
472             
473              ALLOCATE(tabvar3d(SIZE(tabinterp4d,1),SIZE(tabinterp4d,2),SIZE(tabinterp4d,3)))
474              tabvar3d=tabinterp4d(:,:,:,1)
475              CALL Write_Ncdf_var(TRIM(varname),dimnames,Child_file,tabvar3d,1,'double')
476              DEALLOCATE(tabvar3d)
477           ELSE
478              dimnames(3)=vert_coord_name
479              dimnames(4)=TRIM(timedimname)
480             
481              CALL Write_Ncdf_var(TRIM(varname),dimnames,Child_file,tabinterp4d,1,n,'double')
482           ENDIF
483           !             
484           !
485           CALL Copy_Ncdf_att(TRIM(varname),TRIM(restart_file),Child_file)
486           !
487           !
488           IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,src_add,dst_add)   
489           !
490        END DO
491        !
492        DEALLOCATE(detected_pts)
493        DEALLOCATE(tabinterp4d)
494        DEALLOCATE(tabvar1,tabvar2,tabvar3)
495        DEALLOCATE(masksrc) 
496        !                     
497     ENDIF
498
499  END DO
500
501  ! change the before fields
502  IF(rhot == 1) THEN
503     WRITE(*,*) ''   
504     WRITE(*,*) 'no time interpolation (time refinement ratio = 1)'
505  ELSE   
506     now_wght = (rhot-1.)/rhot
507     before_wght = 1./rhot
508     !   
509     ! --- 4D variables --- !               
510     CALL Read_Ncdf_var('un',Child_file,un)     
511     CALL Read_Ncdf_var('vn',Child_file,vn)
512     CALL Read_Ncdf_var('ub',Child_file,ub)
513     CALL Read_Ncdf_var('vb',Child_file,vb)
514     ub = now_wght*un + before_wght*ub
515     vb = now_wght*vn + before_wght*vb
516     !
517     CALL Read_Ncdf_var('tb',Child_file,tb)
518     CALL Read_Ncdf_var('tn',Child_file,tn)
519     tb = now_wght*tn + before_wght*tb
520     !
521     CALL Read_Ncdf_var('sb',Child_file,sb)
522     CALL Read_Ncdf_var('sn',Child_file,sn)           
523     sb = now_wght*sn + before_wght*sb
524     !
525     CALL Read_Ncdf_var('e3t_b',Child_file,e3t_b)
526     CALL Read_Ncdf_var('e3t_n',Child_file,e3t_n)           
527     e3t_b = now_wght*e3t_n + before_wght*e3t_b
528     !
529     dimnames(1)='x'
530     dimnames(2)='y'
531     dimnames(3)='nav_lev'
532     dimnames(4)=TRIM(timedimname)
533     CALL Write_Ncdf_var('ub',dimnames,Child_file,ub,'double')
534     CALL Write_Ncdf_var('vb',dimnames,Child_file,vb,'double')
535     CALL Write_Ncdf_var('tb',dimnames,Child_file,tb,'double')
536     CALL Write_Ncdf_var('sb',dimnames,Child_file,sb,'double')
537     CALL Write_Ncdf_var('e3t_b',dimnames,Child_file,e3t_b,'double')
538     !
539     DEALLOCATE(un,ub,vn,vb,tn,tb,sn,sb,e3t_n,e3t_b)
540     !----------------------               
541     !
542     ! --- 3D variables --- !       
543     CALL Read_Ncdf_var('sshb',Child_file,sshb)
544     CALL Read_Ncdf_var('sshn',Child_file,sshn)
545     sshb = now_wght*sshn + before_wght*sshb
546     !
547     dimnames(1)='x'
548     dimnames(2)='y'
549     dimnames(3)=TRIM(timedimname)
550     CALL Write_Ncdf_var('sshb',dimnames,Child_file,sshb,'double')
551     !               
552     DEALLOCATE(sshn,sshb) 
553     !----------------------               
554     !   
555  ENDIF
556  !
557  !
558  WRITE(*,*) ' '
559  WRITE(*,*) ' --- list of all the variables that have been interpolated --- '
560  WRITE(*,*) Ncdf_varname
561  WRITE(*,*) ' '
562  WRITE(*,*) '******* restart file successfully created *******' 
563  WRITE(*,*) ' '
564  !
565  STOP
566END PROGRAM
567
568
Note: See TracBrowser for help on using the repository browser.