source: utils/tools/NESTING/src/agrif_create_restart_trc.f90 @ 10398

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

BGC restart should now work in the nesting tool. All of that must of course be tested in realistic simulations. Waiting for the feedbacks…

  • Property svn:keywords set to Id
File size: 13.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_trc
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_TRC                   *
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*3 :: prefix
42  CHARACTER*20:: suffix 
43  CHARACTER*100 :: Child_file,Childcoordinates,varname,varname2,Child_Bathy_Level,Child_Bathy_Meter   
44  REAL*8, POINTER, DIMENSION(:,:,:) :: tabvar3d => NULL()
45  REAL*8, POINTER, DIMENSION(:,:,:,:) :: trb,trn => NULL()
46  REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabinterp4d,tabvar1,tabvar2,tabvar3 => NULL()
47  REAL*8, POINTER, DIMENSION(:) :: tabtemp1D,nav_lev => NULL()
48  REAL*8, POINTER, DIMENSION(:,:) :: tabtemp2D => NULL()
49  INTEGER,DIMENSION(:),POINTER :: src_add,dst_add  => NULL()
50  REAL*8,DIMENSION(:,:),POINTER :: matrix => NULL()
51  LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL()
52  LOGICAL, DIMENSION(:,:,:), POINTER :: detected_pts
53  LOGICAL :: Interpolation,Extrapolation,Pacifique
54  INTEGER :: narg,iargc,ncid,x,y,z,nbvert_lev
55  REAL*8 :: now_wght,before_wght
56  INTEGER :: ii,ji,jj,jk,status,varid,numdims
57  CHARACTER(len=20),DIMENSION(4) :: dimnames
58  CHARACTER(len=80) :: namelistname
59  TYPE(Coordinates) :: G0,G1
60  REAL*8 :: tabtemp0dreal
61
62  LOGICAL, PARAMETER :: conservation = .FALSE.
63  !       
64  narg = iargc()
65  IF (narg == 0) THEN
66     namelistname = 'namelist.input'
67  ELSE
68     CALL getarg(1,namelistname)
69  ENDIF
70  !
71  ! read input file
72  !
73  CALL read_namelist(namelistname)
74  !     
75  IF(TRIM(restart_trc_file) == '') THEN
76     WRITE(*,*) 'no tracers restart file specified in ',TRIM(namelistname)
77     STOP
78  ENDIF
79
80  !
81  WRITE(*,*) ''
82  WRITE(*,*) 'Interpolation of restart file : ',TRIM(restart_trc_file)
83  WRITE(*,*) ''
84  !
85  CALL Read_Ncdf_VarName(restart_trc_file,Ncdf_varname)
86  !       
87  CALL set_child_name(parent_coordinate_file,Childcoordinates)   
88  IF( TRIM(parent_bathy_level) /= '' )   CALL set_child_name(parent_bathy_level,Child_Bathy_Level) 
89  IF( TRIM(parent_bathy_meter) /= '' )   CALL set_child_name(parent_bathy_meter,Child_Bathy_Meter)   
90  !
91  ! create this file
92  !
93  CALL set_child_name(restart_trc_file,Child_file)
94  status = nf90_create(Child_file,NF90_WRITE,ncid)
95  status = nf90_close(ncid)
96  WRITE(*,*) 'Child grid restart file name = ',TRIM(Child_file)     
97  WRITE(*,*) ''
98  !
99  ! read dimensions in parent restart file
100  !
101  CALL Read_Ncdf_dim('x',restart_trc_file,x)
102  CALL Read_Ncdf_dim('y',restart_trc_file,y) 
103  CALL Read_Ncdf_dim('nav_lev',restart_trc_file,z)
104
105  IF( z .NE. N ) THEN
106     WRITE(*,*) '***'
107     WRITE(*,*) 'Number of vertical levels doesn t match between namelist and restart file'
108     WRITE(*,*) 'Please check the values in namelist file'
109     STOP
110  ENDIF
111  !
112  ! mask initialization for extrapolation and interpolation
113  !
114  WRITE(*,*) 'mask initialisation on coarse and fine grids'
115  !           
116  status = Read_Local_Coordinates(parent_coordinate_file,G0,(/jpizoom,jpjzoom/),(/x,y/))
117  status = Read_Coordinates(Childcoordinates,G1,Pacifique)
118  !
119  !longitude modification if child domain covers Pacific ocean area
120  !     
121  IF( Pacifique ) THEN
122     WHERE( G0%nav_lon < 0 )
123        G0%nav_lon = G0%nav_lon + 360.
124     END WHERE
125     WHERE( G1%nav_lon < 0 )
126        G1%nav_lon = G1%nav_lon + 360.
127     END WHERE
128  ENDIF
129  !
130  ! one needs bathy_level
131  IF( TRIM(parent_bathy_level) /= '' ) THEN
132     status = Read_bathy_level(TRIM(parent_bathy_level),G0)
133     status = Read_bathy_level(TRIM(child_bathy_level),G1)
134  ELSE
135     status = read_bathy_meter(TRIM(parent_bathy_meter),G0)
136     status = read_bathy_meter(TRIM(child_bathy_meter),G1)
137     CALL meter_to_levels(G0)
138     CALL meter_to_levels(G1)
139  ENDIF
140  ! get masks
141  CALL Init_tmask(parent_bathy_level,G0,x,y)
142  CALL Init_tmask(child_bathy_level,G1,nxfin,nyfin)
143 
144!!clem  G0%tmask = 1.   
145
146  ! which dataset
147  status = nf90_open(TRIM(restart_trc_file), NF90_NOWRITE, ncid) ! Open dataset
148  DO jk = 1, z
149     ALLOCATE(tabvar1(x,y,1,1))
150     !
151     status = nf90_inq_varid(ncid, "TRNDIC", VarId) !PISCES
152     IF (status == nf90_noerr) THEN
153        CALL Read_Ncdf_var('TRNDIC',TRIM(restart_trc_file),tabvar1,1,jk)
154     ELSE
155        status = nf90_inq_varid(ncid, "TRNNO3"  , VarId) ! LOBSTER
156        IF (status == nf90_noerr) THEN
157           CALL Read_Ncdf_var('TRNNO3',TRIM(restart_trc_file),tabvar1,1,jk)
158        ELSE
159           status = nf90_inq_varid(ncid, "TRNCFC11", VarId) ! CFC
160           IF (status == nf90_noerr) THEN
161              CALL Read_Ncdf_var('TRNCFC11',TRIM(restart_trc_file),tabvar1,1,jk)
162           ELSE
163              status = nf90_inq_varid(ncid, "TRNCLR  ", VarId) ! My TRC
164              IF (status == nf90_noerr) THEN
165                 CALL Read_Ncdf_var('TRNCLR',TRIM(restart_trc_file),tabvar1,1,jk)
166              ELSE
167                 WRITE(*,*) 'No suitable tracer found to build the mask '
168              ENDIF
169           ENDIF
170        ENDIF
171     ENDIF
172     WHERE( tabvar1(:,:,1,1) == 0. ) 
173        G0%tmask(:,:,jk) = 0.
174     END WHERE
175     DEALLOCATE(tabvar1)
176  END DO
177
178  !
179  ! write dimensions in output file
180  WRITE(*,*) 'write dimensions'
181  !         
182  CALL Write_Ncdf_dim('x',Child_file,nxfin)
183  CALL Write_Ncdf_dim('y',Child_file,nyfin)
184  CALL Write_Ncdf_dim('nav_lev',Child_file,z)
185  CALL Write_Ncdf_dim('time_counter',Child_file,0) 
186  !
187  !
188  DO ii = 1,SIZE(Ncdf_varname)     
189     !     
190     ! loop on variables names
191     varname = TRIM(Ncdf_varname(ii))
192     WRITE(*,*) 'var = ',TRIM(varname)     
193     !     
194     SELECT CASE (TRIM(varname))
195        !
196     CASE('nav_lon')
197        CALL Read_Ncdf_var('nav_lon',TRIM(Childcoordinates),tabtemp2D) 
198        CALL Write_Ncdf_var('nav_lon',(/'x','y'/),Child_file,tabtemp2D,'float')
199        CALL Copy_Ncdf_att('nav_lon',TRIM(restart_trc_file),Child_file,MINVAL(tabtemp2D),MAXVAL(tabtemp2D))
200        DEALLOCATE(tabtemp2D)
201        Interpolation = .FALSE.
202        !       
203     CASE('nav_lat')             
204        CALL Read_Ncdf_var('nav_lat',TRIM(Childcoordinates),tabtemp2D) 
205        CALL Write_Ncdf_var('nav_lat',(/'x','y'/),Child_file,tabtemp2D,'float')
206        CALL Copy_Ncdf_att('nav_lat',TRIM(restart_trc_file),Child_file,MINVAL(tabtemp2D),MAXVAL(tabtemp2D)) 
207        DEALLOCATE(tabtemp2D)
208        Interpolation = .FALSE.
209        !
210     CASE('nav_lev')
211        CALL Read_Ncdf_var('nav_lev',TRIM(restart_trc_file),nav_lev) 
212        CALL Write_Ncdf_var('nav_lev','z',Child_file,nav_lev,'float')
213        CALL Copy_Ncdf_att('nav_lev',TRIM(restart_trc_file),Child_file)     
214        Interpolation = .FALSE.
215        !
216     CASE('time_counter')
217        CALL Read_Ncdf_var('time_counter',TRIM(restart_trc_file),tabtemp1D) 
218        tabtemp1D = tabtemp1D * rhot
219        CALL Write_Ncdf_var('time_counter','time_counter',Child_file,tabtemp1D,'double')
220        CALL Copy_Ncdf_att('time_counter',TRIM(restart_trc_file),Child_file) 
221        DEALLOCATE(tabtemp1D)
222        Interpolation = .FALSE.
223        !
224     CASE('kt','ndastp','adatrj','ntime','rdttrc1') 
225        CALL Read_Ncdf_var(TRIM(varname),TRIM(restart_trc_file),tabtemp0dreal) 
226        SELECT CASE (TRIM(varname))
227        CASE('rdttrc1')
228           tabtemp0dreal = tabtemp0dreal / rhot
229        CASE('kt')
230           tabtemp0dreal = tabtemp0dreal * rhot
231        END SELECT
232        CALL Write_Ncdf_var(TRIM(varname),Child_file,tabtemp0dreal,'double')
233        CALL Copy_Ncdf_att(varname,TRIM(restart_trc_file),Child_file) 
234        Interpolation = .FALSE.
235        !
236     CASE DEFAULT
237        IF( Get_NbDims(TRIM(varname),TRIM(restart_trc_file)) == 4 ) THEN
238           vert_coord_name = 'nav_lev'
239        ELSEIF( Get_NbDims(TRIM(varname),TRIM(restart_trc_file)) == 3 ) THEN
240           vert_coord_name = '1'
241        ENDIF
242        posvar='T'
243        Interpolation = .TRUE.           
244        !     
245     END SELECT
246     !     
247     ! --- start interpolation --- !
248     IF( Interpolation ) THEN
249        !       
250        IF( vert_coord_name == '1' ) THEN
251           nbvert_lev = 1
252        ELSE
253           nbvert_lev = z
254        ENDIF
255
256        ALLOCATE(detected_pts(SIZE(G0%tmask,1),SIZE(G0%tmask,2),nbvert_lev))                             
257        ALLOCATE(tabvar1(x,y,1,2))
258        ALLOCATE(tabvar2(x,y,1,1))
259        ALLOCATE(tabvar3(x,y,1,1))
260        ALLOCATE(masksrc(x,y))
261        ALLOCATE(tabinterp4d(nxfin,nyfin,1,1)) 
262
263        !       
264        DO n = 1,nbvert_lev
265           !
266           WRITE(*,*) 'interpolate/extrapolate for vertical level = ',n   
267           !
268           CALL Read_Ncdf_var(varname,TRIM(restart_trc_file),tabvar1,1,n)
269           IF(n==1) THEN
270              !                           
271           ELSE IF (n==2) THEN
272              tabvar2(:,:,:,1) = tabvar1(:,:,:,2)
273           ELSE
274              tabvar3(:,:,:,1) = tabvar2(:,:,:,1)
275              tabvar2(:,:,:,1) = tabvar1(:,:,:,2)
276           ENDIF
277           !
278           !
279           IF(MAXVAL(G1%tmask(:,:,n)) == 0.) THEN
280              tabinterp4d = 0.0
281              WRITE(*,*) 'only land points on level ',n                     
282           ELSE
283              CALL extrap_detect(G0,G1,detected_pts(:,:,n),n)                                           
284
285              CALL correct_field(detected_pts(:,:,n),tabvar1,tabvar2,tabvar3,G0,nav_lev,masksrc,n)                                 
286             
287              ! for the following variables, you do not want to mask the values
288              IF(  TRIM(varname) == 'Silicalim' .OR. TRIM(varname) == 'Silicamax' ) THEN
289                 masksrc(:,:) = .TRUE.
290              ENDIF
291             
292              SELECT CASE(TRIM(interp_type))
293              CASE('bilinear')                                                       
294                 CALL get_remap_matrix(G0%nav_lat,G1%nav_lat, &
295                      G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
296                 CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
297                      matrix,src_add,dst_add)     
298              CASE('bicubic')                                   
299                 CALL get_remap_bicub(G0%nav_lat,G1%nav_lat, &
300                      G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
301                 CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),&
302                      nxfin,nyfin,matrix,src_add,dst_add) 
303              END SELECT
304              !                     
305              IF( conservation ) THEN ! clem: it currently does not work
306                 CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
307                    G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
308              ENDIF
309           ENDIF
310
311           IF( ALL(masksrc) ) THEN
312              tabinterp4d(:,:,1,1) =  tabinterp4d(:,:,1,1)
313           ELSE
314              tabinterp4d(:,:,1,1) =  tabinterp4d(:,:,1,1) * G1%tmask(:,:,n)
315           ENDIF
316
317           !
318           dimnames(1)='x'
319           dimnames(2)='y'
320           IF( vert_coord_name == '1' ) THEN
321              dimnames(3)='time_counter'
322             
323              ALLOCATE(tabvar3d(SIZE(tabinterp4d,1),SIZE(tabinterp4d,2),SIZE(tabinterp4d,3)))
324              tabvar3d=tabinterp4d(:,:,:,1)
325              CALL Write_Ncdf_var(TRIM(varname),dimnames,Child_file,tabvar3d,1,'double')
326              DEALLOCATE(tabvar3d)
327           ELSE
328              dimnames(3)=vert_coord_name
329              dimnames(4)='time_counter'
330             
331              CALL Write_Ncdf_var(TRIM(varname),dimnames,Child_file,tabinterp4d,1,n,'double')
332           ENDIF
333           !             
334           !
335           CALL Copy_Ncdf_att(TRIM(varname),TRIM(restart_trc_file),Child_file)
336           !
337           !
338           IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,src_add,dst_add)   
339           !
340           !
341        END DO
342        !
343        DEALLOCATE(detected_pts)
344        DEALLOCATE(tabinterp4d)
345        DEALLOCATE(tabvar1,tabvar2,tabvar3)
346        DEALLOCATE(masksrc) 
347        !                     
348     ENDIF
349
350     ! change the before fields
351     prefix = varname(1:3)
352     suffix = varname(4:LEN_TRIM(varname))   
353
354     IF(rhot == 1 .OR. prefix/= 'TRB') THEN
355        WRITE(*,*) ''   
356        WRITE(*,*) 'no time interpolation for ',TRIM(varname)
357     ELSE   
358        ALLOCATE(trn(nxfin,nyfin,z,1),trb(nxfin,nyfin,z,1))
359        varname2 = 'TRN'//TRIM(suffix)
360        now_wght = (rhot-1.)/rhot
361        before_wght = 1./rhot
362        !   
363        CALL Read_Ncdf_var(TRIM(varname),Child_file,trb)
364        CALL Read_Ncdf_var(TRIM(varname2),Child_file,trn)           
365        trb = now_wght*trn + before_wght*trb
366
367        dimnames(1)='x'
368        dimnames(2)='y'
369        dimnames(3)='nav_lev'
370        dimnames(4)='time_counter'
371        CALL Write_Ncdf_var(TRIM(varname),dimnames,Child_file,trb,'double')
372        DEALLOCATE(trn,trb)
373        !   
374     ENDIF
375
376  END DO
377  !
378  WRITE(*,*) ' '
379  WRITE(*,*) ' --- list of all the variables that have been interpolated --- '
380  WRITE(*,*) Ncdf_varname
381  WRITE(*,*) ' '
382  WRITE(*,*) '******* restart file successfully created *******' 
383  WRITE(*,*) ' '
384  !
385  STOP
386END PROGRAM create_rstrt_trc
387
388
Note: See TracBrowser for help on using the repository browser.