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

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

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

Last change on this file since 10384 was 10384, checked in by clem, 5 years ago

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