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

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

attempt to correct several bugs in the NESTING tools. With this version, domcfg.nc file should be written correctly and the partial steps should be taken into account.

  • Property svn:keywords set to Id
File size: 33.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
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  !
20  IMPLICIT NONE
21  !
22  !************************************************************************
23  !                           *
24  ! PROGRAM  CREATE_RSTRT                    *
25  !                           *
26  ! program to interpolate parent grid restart file to child grid    *
27  !                           *
28  !                           *
29  !Interpolation is carried out using bilinear interpolation      *
30  !routine from SCRIP package                *     
31  !                           *
32  !http://climate.lanl.gov/Software/SCRIP/            *
33  !************************************************************************
34  !
35  ! variables declaration
36  !     
37  CHARACTER*20,DIMENSION(:),POINTER :: Ncdf_varname => NULL()
38  CHARACTER*20 :: vert_coord_name
39  CHARACTER*1 :: posvar
40  CHARACTER*100 :: Child_file,Childcoordinates,varname,Childbathy,Childbathymeter   
41  REAL*8, POINTER, DIMENSION(:,:) :: lonChild,latChild => NULL()
42  REAL*8, POINTER, DIMENSION(:,:) :: lonParent,latParent => NULL()
43  REAL*8, POINTER, DIMENSION(:,:,:) :: tabvar3d,tabinterp3d,mask => NULL()
44  REAL*8, POINTER, DIMENSION(:,:,:) :: fmask,fse3u,fse3v,fse3t => NULL()
45  REAL*8, POINTER, DIMENSION(:,:,:,:) :: un,ub,vn,vb,tn,tb => NULL()
46  REAL*8, POINTER, DIMENSION(:,:,:,:) :: sshn,sshb,sb,sn,gcx,gcxb => NULL()     
47  REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabinterp4d,tabvar1,tabvar2,tabvar3 => NULL()
48  REAL*8, POINTER, DIMENSION(:,:,:,:) :: rotn,rotb,hdivn,hdivb => 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,t,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  !     
75  ! Variables for dimg
76  !
77  INTEGER :: ino0, it0, ipcg0, isor0, itke0,nfice,nfbulk
78  INTEGER :: irecl8, irec,ndastp,narea,nsolv
79  INTEGER :: jk,kt            ! dummy loop indices
80  INTEGER :: inum             ! temporary logical unit
81  INTEGER :: ios1 , ios2      ! flag for ice and bulk in the current run
82  INTEGER :: ios3             ! flag for free surface.  0 = none 1 = yes.  0 = none 1 = yes
83  INTEGER :: ios4             ! flag for coupled (1) or not (0)   
84  !       
85  narg = iargc()
86  IF (narg == 0) THEN
87     namelistname = 'namelist.input'
88  ELSE
89     CALL getarg(1,namelistname)
90  ENDIF
91  !
92  ! read input file
93  !
94  CALL read_namelist(namelistname)
95  !     
96  IF(TRIM(restart_file) == '') THEN
97     WRITE(*,*) 'no restart file specified in ',TRIM(namelistname)
98     STOP
99  END IF
100
101  IF (iom_activated) THEN
102     timedimname = 't'
103  ELSE
104     timedimname='time'
105  ENDIF
106
107  !
108  WRITE(*,*) ''
109  WRITE(*,*) 'Interpolation of restart file : ',TRIM(restart_file)
110  WRITE(*,*) ''
111  !
112  CALL Read_Ncdf_VarName(restart_file,Ncdf_varname)
113  !       
114  CALL set_child_name(parent_coordinate_file,Childcoordinates)   
115  CALL set_child_name(parent_bathy_level,Childbathy) 
116  CALL set_child_name(parent_bathy_meter,Childbathymeter)   
117  !
118  ! create this file
119  !
120  IF( .NOT. dimg ) THEN           
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  ENDIF
127
128  !
129  ! read dimensions in parent restart file
130  !
131  CALL Read_Ncdf_dim('x',restart_file,x)
132  CALL Read_Ncdf_dim('y',restart_file,y) 
133  CALL Read_Ncdf_dim('z',restart_file,z)
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
139  IF( z .NE. N ) THEN
140     WRITE(*,*) '***'
141     WRITE(*,*) 'Number of vertical levels doesn t match between namelist and restart file'
142     WRITE(*,*) 'Please check the values in namelist file'
143     STOP
144  ENDIF
145  !
146  ! mask initialization for extrapolation and interpolation
147  !
148  WRITE(*,*) 'mask initialisation on coarse and fine grids'
149  !           
150  status = Read_Local_Coordinates(parent_coordinate_file,G0,(/jpizoom,jpjzoom/),(/x,y/))
151  status = Read_Coordinates(Childcoordinates,G1,Pacifique)
152  !
153  !longitude modification if child domain covers Pacific ocean area
154  !     
155  IF( Pacifique ) THEN
156     !
157     WHERE( G0%nav_lon < 0 )
158        G0%nav_lon = G0%nav_lon + 360.
159     END WHERE
160     !             
161     WHERE( G1%nav_lon < 0 )
162        G1%nav_lon = G1%nav_lon + 360.
163     END WHERE
164     !
165  ENDIF
166  !
167  CALL Init_mask(parent_bathy_level,G0,x,y)
168  CALL Init_mask(childbathy,G1,1,1)
169
170  G0%tmask = 1.   
171
172  DO k=1,z
173     ALLOCATE(tabvar1(x,y,1,1))
174     CALL Read_Ncdf_var('sn',TRIM(restart_file),tabvar1,1,k)
175     WHERE( tabvar1(:,:,1,1) == 0. ) 
176        G0%tmask(:,:,k) = 0.
177     END WHERE
178     DEALLOCATE(tabvar1)
179  END DO
180  !
181  G0%umask(1:x-1,:,:) = G0%tmask(1:x-1,:,:)*G0%tmask(2:x,:,:)
182  G0%vmask(:,1:y-1,:) = G0%tmask(:,1:y-1,:)*G0%tmask(:,2:y,:)
183  !
184  G0%umask(x,:,:) = G0%tmask(x,:,:)
185  G0%vmask(:,y,:) = G0%tmask(:,y,:)
186  !     
187  G0%fmask(1:x-1,1:y-1,:) = G0%tmask(1:x-1,1:y-1,:)*G0%tmask(2:x,1:y-1,:)* &
188       G0%tmask(1:x-1,2:y,:)*G0%tmask(2:x,2:y,:) 
189  !
190  G0%fmask(x,:,:) = G0%tmask(x,:,:)
191  G0%fmask(:,y,:) = G0%tmask(:,y,:)
192  !
193  !   *****   ***  **   ** ******
194  !   *    *   *   * * * * *
195  !   *    *   *   *  *  * *  ***
196  !   *    *   *   *     * *    *
197  !   *****   ***  *     * ******
198  !
199  IF(dimg) THEN
200     !       
201     WRITE(*,*) 'create dimg restart file'
202     DO m = 7,1000 
203        INQUIRE(Unit=inum,Opened=op)
204        IF( .NOT. op ) THEN
205           inum = m
206           EXIT
207        ENDIF
208     ENDDO
209     !
210     inum = 11
211     irecl8 = nxfin * nyfin * 8
212     OPEN(inum,FILE=TRIM(dimg_output_file),FORM='UNFORMATTED',  &
213          ACCESS='DIRECT',RECL=irecl8)
214     !       
215     CALL Read_Ncdf_var('info',TRIM(restart_file),tabtemp4D)   
216     !       
217     ino0 = NINT(tabtemp4D(1,1,1,1))
218     it0 = NINT(tabtemp4D(1,1,2,1))*rhot
219     WRITE(*,*) 'restart file created for kt = ',it0
220     ipcg0 = NINT(tabtemp4D(1,1,3,1))
221     isor0 = NINT(tabtemp4D(1,1,4,1))
222     itke0 = 0
223     IF(tabtemp4D(1,1,5,1)==1.) itke0 = 1 
224     ndastp = NINT(tabtemp4D(1,1,6,1))
225     ! number of elapsed days since the begining of the run             
226     !
227     DEALLOCATE(tabtemp4D)
228     !
229     IF (isor0 + 1 == 3) THEN
230        isor0 = 2
231        ipcg0 = 2
232     ENDIF
233     !       
234     CALL Read_Ncdf_var('nfice',TRIM(restart_file),tabtemp4D)     
235     nfice = NINT(tabtemp4D(1,1,1,1))
236     DEALLOCATE(tabtemp4D)
237     !
238     CALL Read_Ncdf_var('nfbulk',TRIM(restart_file),tabtemp4D)     
239     nfbulk = NINT(tabtemp4D(1,1,1,1))
240     DEALLOCATE(tabtemp4D)
241     !
242     nfice = 5
243     nfbulk = 5
244     !
245     ios1 = 0
246     ios2 = 0
247     ios3 = 1          ! flag for free surface.  0 = none 1 = yes.  0
248     ios4 = 0
249     narea = 1
250     jpi = nxfin
251     jpj = nyfin
252     jpk = z
253     jpni = 1
254     jpnj = 1
255     jpnij = 1
256     narea = 1
257     jpiglo = nxfin
258     jpjglo = nyfin
259     nlcit = nxfin
260     nlcjt = nyfin
261     nldit = 1
262     nldjt = 1
263     nleit = nxfin
264     nlejt = nyfin
265     nimppt = 1
266     njmppt = 1
267     !
268     PRINT*,'jpi = ',jpi
269     PRINT*,'jpj = ',jpj
270     PRINT*,'jpk = ',jpk
271     PRINT*,'nfice = ',nfice
272     PRINT*,'nfbulk = ',nfbulk
273     PRINT*,'ndastp = ',ndastp
274     !
275     WRITE(inum,REC=1) irecl8, ino0, it0, isor0, ipcg0, itke0, &
276          nfice, nfbulk , ios1, ios2, ios3, ios4, &
277          ndastp, adatrj, jpi, jpj, jpk,  &
278          jpni, jpnj, jpnij, narea, jpiglo, jpjglo, &
279          nlcit, nlcjt, nldit, nldjt, nleit, nlejt, nimppt, njmppt
280     !         
281     irec = 2
282     !
283     !   *****    *****    ******
284     !   *        *    *   *   
285     !   *        *    *   ***
286     !   *        *    *   *   
287     !   *****    *****    * 
288     !
289  ELSE
290
291
292     !
293     ! write dimensions in output file
294     !         
295     CALL Write_Ncdf_dim('x',Child_file,nxfin)
296     CALL Write_Ncdf_dim('y',Child_file,nyfin)
297     CALL Write_Ncdf_dim('z',Child_file,z)
298     CALL Write_Ncdf_dim(TRIM(timedimname),Child_file,0) 
299     IF (.NOT.iom_activated) THEN
300        CALL Write_Ncdf_dim('x_a',Child_file,x_a)
301        CALL Write_Ncdf_dim('y_a',Child_file,y_a)
302        CALL Write_Ncdf_dim('z_a',Child_file,z_a)
303        CALL Write_Ncdf_dim('z_b',Child_file,z_b)
304     ENDIF
305     !
306     !
307  ENDIF
308
309
310
311
312
313
314  !
315  !
316  !
317  DO i = 1,SIZE(Ncdf_varname)     
318     !     
319     ! loop on variables names
320     !     
321     SELECT CASE (TRIM(Ncdf_varname(i)))
322        !
323        !copy nav_lon from child coordinates to output file     
324        !
325     CASE('nav_lon')
326        IF(.NOT. dimg ) THEN
327           WRITE(*,*) 'copy nav_lon'
328           CALL Read_Ncdf_var('nav_lon',TRIM(Childcoordinates),tabtemp2D) 
329           CALL Write_Ncdf_var('nav_lon',(/'x','y'/),Child_file,tabtemp2D,'float')
330           CALL Copy_Ncdf_att('nav_lon',TRIM(restart_file),Child_file, &
331                MINVAL(tabtemp2D),MAXVAL(tabtemp2D))
332           DEALLOCATE(tabtemp2D)
333           Interpolation = .FALSE.
334        ENDIF
335        !       
336        !copy nav_lat from child coordinates to output file
337        !
338     CASE('nav_lat')             
339        IF(.NOT. dimg ) THEN
340           WRITE(*,*) 'copy nav_lat'
341           CALL Read_Ncdf_var('nav_lat',TRIM(Childcoordinates),tabtemp2D) 
342           CALL Write_Ncdf_var('nav_lat',(/'x','y'/),Child_file,tabtemp2D,'float')
343           CALL Copy_Ncdf_att('nav_lat',TRIM(restart_file),Child_file, &
344                MINVAL(tabtemp2D),MAXVAL(tabtemp2D)) 
345           DEALLOCATE(tabtemp2D)
346           Interpolation = .FALSE.
347        ENDIF
348        !
349        !copy nav_lev from restart_file to output file
350        !
351     CASE('nav_lev')
352
353        WRITE(*,*) 'copy nav_lev'
354        CALL Read_Ncdf_var('nav_lev',TRIM(restart_file),nav_lev) 
355        IF(.NOT. dimg ) THEN
356           CALL Write_Ncdf_var('nav_lev','z',Child_file,nav_lev,'float')
357           CALL Copy_Ncdf_att('nav_lev',TRIM(restart_file),Child_file)     
358        ENDIF
359        Interpolation = .FALSE.
360        !
361        !copy time from restart_file to output file                       
362        !
363     CASE('time')
364        IF(.NOT. dimg ) THEN
365           WRITE(*,*) 'copy time'
366           CALL Read_Ncdf_var('time',TRIM(restart_file),tabtemp1D) 
367           CALL Write_Ncdf_var('time',TRIM(timedimname),Child_file,tabtemp1D,'float')
368           CALL Copy_Ncdf_att('time',TRIM(restart_file),Child_file) 
369           DEALLOCATE(tabtemp1D)
370           Interpolation = .FALSE.
371        ENDIF
372        !copy time from restart_file to output file                       
373        !
374     CASE('time_counter')
375        IF(.NOT. dimg ) THEN
376           WRITE(*,*) 'copy time_counter'
377           CALL Read_Ncdf_var('time_counter',TRIM(restart_file),tabtemp1D) 
378           tabtemp1D = tabtemp1D * rhot
379           CALL Write_Ncdf_var('time_counter',TRIM(timedimname),Child_file,tabtemp1D,'double')
380           CALL Copy_Ncdf_att('time_counter',TRIM(restart_file),Child_file) 
381           DEALLOCATE(tabtemp1D)
382           Interpolation = .FALSE.
383        ENDIF
384        !
385        !copy time_steps from restart_file to output file             
386        !
387     CASE('time_steps')
388        IF(.NOT. dimg ) THEN
389           WRITE(*,*) 'copy time_steps'
390           CALL Read_Ncdf_var('time_steps',TRIM(restart_file),tabtemp1DInt) 
391           CALL Write_Ncdf_var('time_steps',TRIM(timedimname),Child_file,tabtemp1DInt,'integer')
392           CALL Copy_Ncdf_att('time_steps',TRIM(restart_file),Child_file) 
393           DEALLOCATE(tabtemp1DInt)
394           Interpolation = .FALSE.   
395        ENDIF
396        !
397        !copy info from restart_file to output file
398        !
399     CASE('info') 
400        IF(.NOT. dimg ) THEN
401           WRITE(*,*) 'copy info'       
402           CALL Read_Ncdf_var('info',TRIM(restart_file),tabtemp4D) 
403           dimnames(1)='x_a'
404           dimnames(2)='y_a'
405           dimnames(3)='z_a'
406           dimnames(4)=TRIM(timedimname)           
407           CALL Write_Ncdf_var('info',dimnames,Child_file,tabtemp4D,'double')
408           CALL Copy_Ncdf_att('info',TRIM(restart_file),Child_file) 
409           DEALLOCATE(tabtemp4D)
410           Interpolation = .FALSE. 
411        ENDIF
412        !
413     CASE('nfice','nfbulk','kt','ndastp','adatrj','rdt') 
414        IF(.NOT. dimg ) THEN
415           WRITE(*,*) 'copy ',TRIM(Ncdf_varname(i))
416           IF (iom_activated) THEN
417              CALL Read_Ncdf_var(TRIM(Ncdf_varname(i)),TRIM(restart_file),tabtemp0dreal) 
418              SELECT CASE (TRIM(Ncdf_varname(i)))
419              CASE('rdt')
420                 tabtemp0dreal = tabtemp0dreal/rhot
421              CASE('kt')
422                 tabtemp0dreal = tabtemp0dreal * rhot
423              END SELECT
424              CALL Write_Ncdf_var(TRIM(Ncdf_varname(i)),Child_file,tabtemp0dreal,'double')
425           ELSE
426              CALL Read_Ncdf_var(TRIM(Ncdf_varname(i)),TRIM(restart_file),tabtemp4D) 
427              dimnames(1)='x_a'
428              dimnames(2)='y_a'
429              dimnames(3)='z_b'
430              dimnames(4)=TRIM(timedimname)           
431              CALL Write_Ncdf_var(TRIM(Ncdf_varname(i)),dimnames,Child_file,tabtemp4D,'double')
432         DEALLOCATE(tabtemp4D)
433           ENDIF
434           CALL Copy_Ncdf_att(TRIM(Ncdf_varname(i)),TRIM(restart_file),Child_file) 
435           Interpolation = .FALSE.
436        ENDIF
437        !
438        ! Variable interpolation according to their position on grid
439        !                                   
440     CASE('un','ub') 
441        varname = Ncdf_varname(i)
442
443        IF(TRIM(varname)=='un') irec = 6 * z + 2
444        IF(TRIM(varname)=='ub') irec = 2
445
446        WRITE(*,*) TRIM(varname),'interpolation ...'   
447        vert_coord_name = 'z'             
448        posvar='U'
449        Interpolation = .TRUE. 
450        !
451     CASE('u_io') 
452        varname = Ncdf_varname(i) 
453        WRITE(*,*) TRIM(varname),'interpolation ...'   
454        vert_coord_name = 'z_b'             
455        posvar='U'
456        Interpolation = .TRUE. 
457        !                         
458     CASE('vn','vb')
459        varname = Ncdf_varname(i)
460
461        IF(TRIM(varname)=='vn') irec = 7 * z + 2
462        IF(TRIM(varname)=='vb') irec = 2 + z
463
464        WRITE(*,*) TRIM(varname),'interpolation ...'     
465        vert_coord_name = 'z'
466        posvar='V'
467        Interpolation = .TRUE.     
468        !             
469     CASE('v_io')
470        varname = Ncdf_varname(i)
471        WRITE(*,*) TRIM(varname),'interpolation ...'     
472        vert_coord_name = 'z_b'
473        posvar='V'
474        Interpolation = .TRUE.
475        !             
476     CASE('gcx','gcxb','sshb','sshn','sst_io','sss_io','gsst')
477        varname = Ncdf_varname(i)   
478
479        IF(TRIM(varname)=='gcx') irec = 12 * z + 2
480        IF(TRIM(varname)=='gcxb')  irec = 12 * z + 3
481        IF(TRIM(varname)=='sshb') irec = 12 * z + 4
482        IF(TRIM(varname)=='sshn') irec = 12 * z + 5
483
484        WRITE(*,*) TRIM(varname),'interpolation ...' 
485        vert_coord_name = 'z_b'             
486        posvar='T'
487        Interpolation = .TRUE.     
488
489       
490     CASE ('tb','sb','sn','tn')
491        varname = Ncdf_varname(i) 
492
493        IF(TRIM(varname)=='sn') irec = 9 * z + 2
494        IF(TRIM(varname)=='tn') irec = 8 * z + 2
495        IF(TRIM(varname)=='sb') irec = 3 * z + 2
496        IF(TRIM(varname)=='tb') irec = 2 * z + 2
497
498        WRITE(*,*) TRIM(varname),'interpolation ...'     
499        vert_coord_name = 'z'
500        posvar='T'
501        Interpolation = .TRUE.           
502
503     CASE('en')
504        varname = Ncdf_varname(i) 
505        irec = 12 * z + 6
506        WRITE(*,*) TRIM(varname),'interpolation ...'     
507        vert_coord_name = 'z'
508        posvar='T'
509        Interpolation = .TRUE.             
510
511     CASE ('rotb','rotn','hdivb','hdivn')
512        Interpolation = .FALSE.
513        !
514     END SELECT
515     !     
516     IF( Interpolation ) THEN
517        !       
518        IF( vert_coord_name == 'z') THEN
519           nbvert_lev = z 
520        ELSE IF( vert_coord_name == 'z_b') THEN
521           IF (iom_activated) THEN
522              nbvert_lev=1
523           ELSE
524              nbvert_lev = z_b
525           ENDIF
526        END IF
527        !                 
528
529        !
530        t = 1
531
532        ALLOCATE(detected_pts(SIZE(G0%tmask,1),SIZE(G0%tmask,2),nbvert_lev))                             
533        ALLOCATE(tabvar1(x,y,1,2))
534        ALLOCATE(tabvar2(x,y,1,1))
535        ALLOCATE(tabvar3(x,y,1,1))
536        ALLOCATE(masksrc(x,y))
537        ALLOCATE(tabinterp4d(nxfin,nyfin,1,1)) 
538
539        !       
540        DO n = 1,nbvert_lev
541           !
542           WRITE(*,*) 'interpolate/extrapolate ', &
543                TRIM(varname),' for vertical level = ',n   
544           !                           
545           !                                                       
546           !                            If(n==1) then
547           !                                      Call Read_Ncdf_var(varname,TRIM(restart_file),tabvar1,t,n)   
548           !                            else if (n==2) then
549           !                                      Call Read_Ncdf_var(varname,TRIM(restart_file),tabvar2,t,n-1)
550           !                                      Call Read_Ncdf_var(varname,TRIM(restart_file),tabvar1,t,n)           
551           !                            else
552           !                                      Call Read_Ncdf_var(varname,TRIM(restart_file),tabvar3,t,n-2)
553           !                                      Call Read_Ncdf_var(varname,TRIM(restart_file),tabvar2,t,n-1)
554           !                                      Call Read_Ncdf_var(varname,TRIM(restart_file),tabvar1,t,n)
555           !                            endif                                                                           
556           !
557           CALL Read_Ncdf_var(varname,TRIM(restart_file),tabvar1,t,n)
558           IF(n==1) THEN
559              !                           
560           ELSE IF (n==2) THEN
561              tabvar2(:,:,:,1) = tabvar1(:,:,:,2)
562           ELSE
563              tabvar3(:,:,:,1) = tabvar2(:,:,:,1)
564              tabvar2(:,:,:,1) = tabvar1(:,:,:,2)
565
566           ENDIF
567           !                           
568           SELECT CASE(posvar)
569              !
570           CASE('T')
571              !
572              IF(MAXVAL(G1%tmask(:,:,n)) == 0.) THEN           
573                 tabinterp4d = 0.0
574                 WRITE(*,*) 'only land points on level ',n                     
575              ELSE                       
576
577                 CALL extrap_detect(G0,G1,detected_pts(:,:,n),n)                                           
578
579                 CALL correct_field(detected_pts(:,:,n),tabvar1,tabvar2,&
580                      tabvar3,G0,nav_lev,masksrc,n)                                 
581
582                 SELECT CASE(TRIM(interp_type))
583                 CASE('bilinear')                                                       
584                    CALL get_remap_matrix(G0%nav_lat,G1%nav_lat, &
585                         G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
586                    CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
587                         matrix,src_add,dst_add)     
588                 CASE('bicubic')                                   
589                    CALL get_remap_bicub(G0%nav_lat,G1%nav_lat, &
590                         G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
591                    CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),&
592                         nxfin,nyfin,matrix,src_add,dst_add) 
593                 END SELECT
594                 !                     
595                 CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
596                      G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
597
598
599              ENDIF
600
601              tabinterp4d(:,:,1,1) =  tabinterp4d(:,:,1,1) * G1%tmask(:,:,n)
602              !
603           CASE('U')
604              !
605              IF(MAXVAL(G1%umask(:,:,n)) == 0) THEN
606                 tabinterp4d = 0.0
607                 WRITE(*,*) 'only land points on level ',n 
608              ELSE     
609                 !
610                 CALL extrap_detect(G0,G1,detected_pts(:,:,n),n,'U') 
611                 CALL correct_field(detected_pts(:,:,n),tabvar1,tabvar2,&
612                      tabvar3,G0,nav_lev,masksrc,n,'U')
613                 !                                                           
614                 SELECT CASE(TRIM(interp_type))
615                 CASE('bilinear')                                                       
616                    CALL get_remap_matrix(G0%gphiu,G1%gphiu,   &
617                         G0%glamu,G1%glamu,masksrc,matrix,src_add,dst_add)
618                    CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
619                         matrix,src_add,dst_add)     
620                 CASE('bicubic')                                   
621                    CALL get_remap_bicub(G0%gphiu,G1%gphiu,   &
622                         G0%glamu,G1%glamu,masksrc,matrix,src_add,dst_add)
623                    CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),&
624                         nxfin,nyfin,matrix,src_add,dst_add)                       
625                 END SELECT
626                 !                     
627                 CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
628                      G0%e1u,G0%e2u,G1%e1u,G1%e2u,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
629
630              ENDIF
631
632              tabinterp4d(:,:,1,1) =  tabinterp4d(:,:,1,1) * G1%umask(:,:,n)
633              !
634           CASE('V')
635              !     
636              IF(MAXVAL(G1%vmask(:,:,n)) == 0) THEN
637                 tabinterp4d = 0.0
638                 WRITE(*,*) 'only land points on level ',n 
639              ELSE     
640                 !
641
642                 CALL extrap_detect(G0,G1,detected_pts(:,:,n),n,'V')
643
644                 CALL correct_field(detected_pts(:,:,n),tabvar1,tabvar2,&
645                      tabvar3,G0,nav_lev,masksrc,n,'V')
646                 !                                                           
647                 SELECT CASE(TRIM(interp_type))
648                 CASE('bilinear')                                                       
649                    CALL get_remap_matrix(G0%gphiv,G1%gphiv,   &
650                         G0%glamv,G1%glamv,masksrc,matrix,src_add,dst_add)
651                    CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
652                         matrix,src_add,dst_add)     
653                 CASE('bicubic')                                   
654                    CALL get_remap_bicub(G0%gphiv,G1%gphiv,   &
655                         G0%glamv,G1%glamv,masksrc,matrix,src_add,dst_add)
656                    CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),&
657                         nxfin,nyfin,matrix,src_add,dst_add)                       
658                 END SELECT
659                 !                     
660                 CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
661                      G0%e1v,G0%e2v,G1%e1v,G1%e2v,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
662
663              ENDIF
664
665              tabinterp4d(:,:,1,1) =  tabinterp4d(:,:,1,1) * G1%vmask(:,:,n)
666              !
667           END SELECT
668           !
669           IF(dimg) THEN
670              !
671              WRITE(inum,REC=irec) tabinterp4d(:,:,1,1)
672              irec = irec + 1
673              !
674           ELSE
675              !
676              dimnames(1)='x'
677              dimnames(2)='y'
678              IF ((iom_activated).AND.(vert_coord_name == 'z_b')) THEN
679                 dimnames(3)=TRIM(timedimname)
680              ELSE
681                 dimnames(3)=vert_coord_name
682                 dimnames(4)=TRIM(timedimname)
683              ENDIF
684              !             
685              IF(TRIM(varname)=='gcx' .OR. TRIM(varname)=='gcxb') THEN
686                 PRINT*,TRIM(varname),MAXVAL(tabinterp4d),MINVAL(tabinterp4d)
687              ENDIF
688              IF ((iom_activated).AND.(vert_coord_name == 'z_b')) THEN
689                 ALLOCATE(tabvar3d(SIZE(tabinterp4d,1),SIZE(tabinterp4d,2),SIZE(tabinterp4d,3)))
690                 tabvar3d=tabinterp4d(:,:,:,1)
691                 CALL Write_Ncdf_var(TRIM(varname),dimnames, &
692                      Child_file,tabvar3d,t,'double')
693                 DEALLOCATE(tabvar3d)
694              ELSE
695                 CALL Write_Ncdf_var(TRIM(varname),dimnames, &
696                      Child_file,tabinterp4d,t,n,'double')
697              ENDIF
698              !
699              CALL Copy_Ncdf_att(TRIM(varname),TRIM(restart_file),Child_file)
700              !
701           ENDIF
702           !
703           IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,src_add,dst_add)   
704           !
705        END DO
706        !
707        DEALLOCATE(detected_pts)
708        DEALLOCATE(tabinterp4d)
709        DEALLOCATE(tabvar1,tabvar2,tabvar3)
710        DEALLOCATE(masksrc) 
711        !                     
712     ENDIF
713
714  END DO
715  !
716  !
717  !
718  !
719  !
720  !
721  !
722  !
723  !
724  ALLOCATE(rotn(nxfin,nyfin,z,1),rotb(nxfin,nyfin,z,1))
725  ALLOCATE(hdivn(nxfin,nyfin,z,1),hdivb(nxfin,nyfin,z,1))
726  !                         
727  !
728  IF(dimg) THEN         
729     ALLOCATE(un(nxfin,nyfin,z,1),ub(nxfin,nyfin,z,1))
730     ALLOCATE(vn(nxfin,nyfin,z,1),vb(nxfin,nyfin,z,1))
731     irec = 6 * z + 2
732     CALL read_dimg_var(inum,irec,un,z) 
733     irec = 2
734     CALL read_dimg_var(inum,irec,ub,z)
735     irec = 7 * z + 2
736     CALL read_dimg_var(inum,irec,vn,z)
737     irec = 2 + z
738     CALL read_dimg_var(inum,irec,vb,z)                     
739  ELSE
740     CALL Read_Ncdf_var('un',Child_file,un)     
741     CALL Read_Ncdf_var('vn',Child_file,vn)
742     !       
743     CALL Read_Ncdf_var('ub',Child_file,ub)
744     CALL Read_Ncdf_var('vb',Child_file,vb)
745  ENDIF
746  !
747
748  IF(rhot == 1) THEN
749     WRITE(*,*) ''   
750     WRITE(*,*) 'no time interpolation (time refinement ratio = 1)'
751  ELSE   
752     WRITE(*,*) ''
753     WRITE(*,*) 'time interpolation'
754     now_wght = (rhot-1.)/rhot
755     before_wght = 1./rhot
756     !   
757     ub = now_wght*un + before_wght*ub
758     vb = now_wght*vn + before_wght*vb
759     !-----------------------                 
760     IF(dimg) THEN
761        ALLOCATE(tn(nxfin,nyfin,z,1),tb(nxfin,nyfin,z,1))
762        irec = 8 * z + 2
763        CALL read_dimg_var(inum,irec,tn,z)
764        irec = 2 * z + 2
765        CALL read_dimg_var(inum,irec,tb,z)               
766     ELSE
767        CALL Read_Ncdf_var('tb',Child_file,tb)
768        CALL Read_Ncdf_var('tn',Child_file,tn)
769     ENDIF
770     !----------------------               
771     tb = now_wght*tn + before_wght*tb
772
773     IF(dimg) THEN
774        irec = 2 * z + 2
775        CALL write_dimg_var(inum,irec,tb,z)               
776     ELSE
777        dimnames(1)='x'
778        dimnames(2)='y'
779        dimnames(3)='z'
780        dimnames(4)=TRIM(timedimname)
781        CALL Write_Ncdf_var('tb',dimnames,Child_file,tb,'double')
782     ENDIF
783     !
784     DEALLOCATE(tn,tb)
785     !----------------------         
786     IF(dimg) THEN
787        ALLOCATE(sn(nxfin,nyfin,z,1),sb(nxfin,nyfin,z,1))
788        irec = 9 * z + 2
789        CALL read_dimg_var(inum,irec,sn,z)
790        irec = 3 * z + 2
791        CALL read_dimg_var(inum,irec,sb,z)             
792     ELSE         
793        CALL Read_Ncdf_var('sb',Child_file,sb)
794        CALL Read_Ncdf_var('sn',Child_file,sn)           
795     ENDIF
796     !----------------------               
797     sb = now_wght*sn + before_wght*sb
798
799     IF(dimg) THEN
800        irec = 3 * z + 2
801        CALL write_dimg_var(inum,irec,sb,z)               
802     ELSE               
803        dimnames(1)='x'
804        dimnames(2)='y'
805        dimnames(3)='z'
806        dimnames(4)=TRIM(timedimname)
807        CALL Write_Ncdf_var('sb',dimnames,Child_file,sb,'double')
808     ENDIF
809     !----------------------                             
810     DEALLOCATE(sn,sb)
811     !         
812     IF(dimg) THEN
813        ALLOCATE(gcx(nxfin,nyfin,1,1),gcxb(nxfin,nyfin,1,1))
814        irec = 12 * z + 2
815        CALL read_dimg_var(inum,irec,gcx,1)
816        irec = 12 * z + 3
817        CALL read_dimg_var(inum,irec,gcxb,1) 
818     ELSE 
819        CALL Read_Ncdf_var('gcx',Child_file,gcx)
820        CALL Read_Ncdf_var('gcxb',Child_file,gcxb)
821     ENDIF
822     !----------------------               
823     gcxb = now_wght*gcx + before_wght*gcxb
824
825     IF(dimg) THEN
826        irec = 12 * z + 3
827        CALL write_dimg_var(inum,irec,gcxb,1)               
828     ELSE               
829        dimnames(1)='x'
830        dimnames(2)='y'
831        dimnames(3)='z_b'
832        dimnames(4)=TRIM(timedimname)
833        CALL Write_Ncdf_var('gcxb',dimnames,Child_file,gcxb,'double')
834     ENDIF
835     !-----------------------               
836     PRINT*,' gcx = ',MAXVAL(gcx),MINVAL(gcx)
837     PRINT*,' gcxb = ',MAXVAL(gcxb),MINVAL(gcxb)
838
839     DEALLOCATE(gcx,gcxb)
840     !
841     IF(dimg) THEN
842        ALLOCATE(sshn(nxfin,nyfin,1,1),sshb(nxfin,nyfin,1,1))
843        irec = 12 * z + 5
844        CALL read_dimg_var(inum,irec,sshn,1)
845        irec = 12 * z + 4
846        CALL read_dimg_var(inum,irec,sshb,1) 
847     ELSE         
848        CALL Read_Ncdf_var('sshb',Child_file,sshb)
849        CALL Read_Ncdf_var('sshn',Child_file,sshn)
850     ENDIF
851     !----------------------               
852     sshb = now_wght*sshn + before_wght*sshb
853
854     IF(dimg) THEN
855        irec = 12 * z + 4
856        CALL write_dimg_var(inum,irec,sshb,1)               
857     ELSE         
858        dimnames(1)='x'
859        dimnames(2)='y'
860        dimnames(3)='z_b'
861        dimnames(4)=TRIM(timedimname)
862        CALL Write_Ncdf_var('sshb',dimnames,Child_file,sshb,'double')
863     ENDIF
864     !               
865     DEALLOCATE(sshb,sshn) 
866
867     !   
868  ENDIF
869  !
870  WRITE(*,*) 'Compute hdivn,rotn with new interpolated fields ...'
871  !
872  !fmask:land/ocean mask at f-point (=0. or 1.)=shlat along lateral boundaries
873  !
874  ALLOCATE(fse3u(nxfin,nyfin,z),fse3v(nxfin,nyfin,z),fse3t(nxfin,nyfin,z)) 
875  !     
876  IF(partial_steps) THEN
877     status = Read_Bathy_Meter(TRIM(Childbathymeter),G1)
878     CALL get_scale_factors( G1,fse3t,fse3u,fse3v )
879  ELSE       
880     fse3t(:,:,:) = 1.0
881     fse3u(:,:,:) = 1.0
882     fse3v(:,:,:) = 1.0
883  ENDIF
884  !
885  !
886  DO k = 1,z 
887     !   
888     DO j = 2,nyfin-1
889        DO i = 2, nxfin-1
890           !
891           hdivn(i,j,k,1) = (G1%e2u(i,j)*fse3u(i,j,k)*un(i,j,k,1)-G1%e2u(i-1,j)*fse3u(i-1,j,k)*un(i-1,j,k,1)      &
892                +    G1%e1v(i,j)*fse3v(i,j,k)*vn(i,j,k,1)-G1%e1v(i,j-1)*fse3v(i,j-1,k)*vn(i,j-1,k,1))   & 
893                / ( G1%e1t(i,j) * G1%e2t(i,j) * fse3t(i,j,k) )
894           !
895           hdivb(i,j,k,1) = (G1%e2u(i,j)*fse3u(i,j,k)*ub(i,j,k,1)-G1%e2u(i-1,j)*fse3u(i-1,j,k)*ub(i-1,j,k,1)      &
896                +    G1%e1v(i,j)*fse3v(i,j,k)*vb(i,j,k,1)-G1%e1v(i,j-1)*fse3v(i,j-1,k)*vb(i,j-1,k,1))   & 
897                / ( G1%e1t(i,j) * G1%e2t(i,j) * fse3t(i,j,k) )
898           !         
899        END DO
900     END DO
901     !
902     !
903     hdivn(1:2,:,k,1) = 0.
904     hdivn(:,1:2,k,1) = 0.
905     hdivn(nxfin-1:nxfin,:,k,1) = 0.
906     hdivn(:,nyfin-1:nyfin,k,1) = 0.
907     hdivb(1:2,:,k,1) = 0.
908     hdivb(:,1:2,k,1) = 0.
909     hdivb(nxfin-1:nxfin,:,k,1) = 0.
910     hdivb(:,nyfin-1:nyfin,k,1) = 0.         
911     !   
912     DO j = 1, nyfin-1
913        DO i = 1, nxfin-1
914           !
915           rotn(i,j,k,1) = (G1%e2v(i+1,j)*vn(i+1,j,k,1)-G1%e2v(i,j)*vn(i,j,k,1)    &
916                - G1%e1u(i,j+1)*un(i,j+1,k,1)+G1%e1u(i,j)*un(i,j,k,1)) &
917                * G1%fmask(i,j,k) / ( G1%e1f(i,j) * G1%e2f(i,j))
918           !
919           rotb(i,j,k,1) = (G1%e2v(i+1,j)*vb(i+1,j,k,1)-G1%e2v(i,j)*vb(i,j,k,1)    &
920                - G1%e1u(i,j+1)*ub(i,j+1,k,1)+G1%e1u(i,j)*ub(i,j,k,1)) &
921                * G1%fmask(i,j,k) / ( G1%e1f(i,j) * G1%e2f(i,j))
922           !             
923        END DO
924     END DO
925     !
926  END DO
927  !
928  PRINT*,' hdivn = ',MAXVAL(hdivn(2:nxfin-1,2:nyfin-1,:,1)),MINVAL(hdivn(2:nxfin-1,2:nyfin-1,:,1))
929  PRINT*,' hdivb = ',MAXVAL(hdivb(2:nxfin-1,2:nyfin-1,:,1)),MINVAL(hdivb(2:nxfin-1,2:nyfin-1,:,1))
930  PRINT*,' rotn  = ',MAXVAL(rotn(2:nxfin-1,2:nyfin-1,:,1)),MINVAL(rotn(2:nxfin-1,2:nyfin-1,:,1))
931  PRINT*,' rotb  = ',MAXVAL(rotb(2:nxfin-1,2:nyfin-1,:,1)),MINVAL(rotb(2:nxfin-1,2:nyfin-1,:,1))
932  !
933  IF(dimg) THEN
934     !
935     irec = 10 * z + 2
936     DO k=1,z
937        WRITE(inum,REC=irec) rotn(:,:,k,1)
938        irec = irec+1
939     ENDDO
940     !         
941     irec = 4 * z + 2
942     DO k=1,z
943        WRITE(inum,REC=irec) rotb(:,:,k,1)
944        irec = irec+1
945     ENDDO
946     !         
947     irec = 11 * z + 2
948     DO k=1,z
949        WRITE(inum,REC=irec) hdivn(:,:,k,1)
950        irec = irec+1
951     ENDDO
952     !         
953     irec = 5 * z + 2
954     DO k=1,z
955        WRITE(inum,REC=irec) hdivb(:,:,k,1)
956        irec = irec+1
957     ENDDO
958     !         
959     CLOSE(inum)
960  ELSE
961     dimnames(1)='x'
962     dimnames(2)='y'
963     dimnames(3)='z'
964     dimnames(4)=TRIM(timedimname)
965     CALL Write_Ncdf_var('rotn',dimnames,Child_file,rotn,'double')
966     CALL Copy_Ncdf_att('rotn',TRIM(restart_file),Child_file)
967     CALL Write_Ncdf_var('hdivn',dimnames,Child_file,hdivn,'double')
968     CALL Copy_Ncdf_att('hdivn',TRIM(restart_file),Child_file)
969     CALL Write_Ncdf_var('rotb',dimnames,Child_file,rotb,'double')
970     CALL Copy_Ncdf_att('rotb',TRIM(restart_file),Child_file)
971     CALL Write_Ncdf_var('hdivb',dimnames,Child_file,hdivb,'double')
972     CALL Copy_Ncdf_att('hdivb',TRIM(restart_file),Child_file)
973  ENDIF
974  !
975  WRITE(*,*) ' '
976  WRITE(*,*) '******* restart file successfully created *******' 
977  WRITE(*,*) ' '
978  !
979  STOP
980END PROGRAM create_rstrt
981
982
Note: See TracBrowser for help on using the repository browser.