source: branches/UKMO/r6232_tracer_advection/NEMOGCM/TOOLS/NESTING/src/agrif_create_restart.f90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

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) == '/NULL') 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_meshmask_file,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_meshmask_file,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)
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','rdttra1') 
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','rdttra1')
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_Bathymeter(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.