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

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

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

  • Property svn:keywords set to Id
File size: 15.1 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  !
20  IMPLICIT NONE
21  !
22  !************************************************************************
23  !                           *
24  ! PROGRAM  CREATE_RSTRT_TRC                   *
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*3 :: prefix
41  CHARACTER*20:: suffix 
42  CHARACTER*100 :: Child_file,Childcoordinates,varname,varname2,Childbathy,Childbathymeter   
43  REAL*8, POINTER, DIMENSION(:,:,:) :: fse3u, fse3v,fse3t => NULL()
44  REAL*8, POINTER, DIMENSION(:,:) :: lonChild,latChild => NULL()
45  REAL*8, POINTER, DIMENSION(:,:) :: lonParent,latParent => NULL()
46  REAL*8, POINTER, DIMENSION(:,:,:) :: cvol => NULL()
47  REAL*8, POINTER, DIMENSION(:,:,:) :: tabvar3d,tabinterp3d,mask => NULL()
48  REAL*8, POINTER, DIMENSION(:,:,:,:) :: trb,trn => NULL()
49  REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabinterp4d,tabvar1,tabvar2,tabvar3 => 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,ztrcor,zvolk
63  INTEGER :: narg,iargc,ncid,x,y,t,z,nbvert_lev,m
64  REAL*8 :: now_wght,before_wght, ztrmas,zcoef
65  INTEGER :: i,j,k,ji,jj,jk,status,varid,numdims
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  !       
76  narg = iargc()
77  IF (narg == 0) THEN
78     namelistname = 'namelist.input'
79  ELSE
80     CALL getarg(1,namelistname)
81  ENDIF
82  !
83  ! read input file
84  !
85  CALL read_namelist(namelistname)
86  !     
87  IF(TRIM(restart_trc_file) == '') THEN
88     WRITE(*,*) 'no tracers restart file specified in ',TRIM(namelistname)
89     STOP
90  ENDIF
91
92  timedimname = 't'
93
94  !
95  WRITE(*,*) ''
96  WRITE(*,*) 'Interpolation of restart file : ',TRIM(restart_trc_file)
97  WRITE(*,*) ''
98  !
99  CALL Read_Ncdf_VarName(restart_trc_file,Ncdf_varname)
100  !       
101  CALL set_child_name(parent_coordinate_file,Childcoordinates)   
102  CALL set_child_name(parent_bathy_level,Childbathy) 
103  CALL set_child_name(parent_bathy_meter,Childbathymeter)   
104  !
105  ! create this file
106  !
107  CALL set_child_name(restart_trc_file,Child_file)
108  status = nf90_create(Child_file,NF90_WRITE,ncid)
109  status = nf90_close(ncid)
110  WRITE(*,*) 'Child grid restart file name = ',TRIM(Child_file)     
111  WRITE(*,*) ''
112  !
113  ! read dimensions in parent restart file
114  !
115  CALL Read_Ncdf_dim('x',restart_trc_file,x)
116  CALL Read_Ncdf_dim('y',restart_trc_file,y) 
117  CALL Read_Ncdf_dim('z',restart_trc_file,z)
118
119  IF( z .NE. N ) THEN
120     WRITE(*,*) '***'
121     WRITE(*,*) 'Number of vertical levels doesn t match between namelist and restart file'
122     WRITE(*,*) 'Please check the values in namelist file'
123     STOP
124  ENDIF
125  !
126  ! mask initialization for extrapolation and interpolation
127  !
128  WRITE(*,*) 'mask initialisation on coarse and fine grids'
129  !           
130  status = Read_Local_Coordinates(parent_coordinate_file,G0,(/jpizoom,jpjzoom/),(/x,y/))
131  status = Read_Coordinates(Childcoordinates,G1,Pacifique)
132  !
133  !longitude modification if child domain covers Pacific ocean area
134  !     
135  IF( Pacifique ) THEN
136     WHERE( G0%nav_lon < 0 )
137        G0%nav_lon = G0%nav_lon + 360.
138     END WHERE
139     WHERE( G1%nav_lon < 0 )
140        G1%nav_lon = G1%nav_lon + 360.
141     END WHERE
142  ENDIF
143  !
144  CALL Init_tmask(parent_bathy_level,G0,x,y)
145  CALL Init_tmask(childbathy,G1,nxfin,nyfin)
146 
147  ALLOCATE(fse3u(nxfin,nyfin,z),fse3v(nxfin,nyfin,z),fse3t(nxfin,nyfin,z),cvol(nxfin,nyfin,z))
148    !     
149    status = Read_Bathy_Meter(TRIM(Childbathymeter),G1)
150    CALL get_scale_factors( G1,fse3t,fse3u,fse3v )
151                         
152  DO jk =1,z
153  cvol(:,:,jk) = G1%e1t(:,:)*G1%e2t(:,:)*fse3t(:,:,jk)*G1%tmask(:,:,jk)
154  END DO 
155!  CALL Init_mask(childbathy,G1,1,1)
156
157  G0%tmask = 1.   
158
159
160  status = nf90_open(TRIM(restart_trc_file), NF90_NOWRITE, ncid) ! Open dataset
161  DO k=1,z
162     ALLOCATE(tabvar1(x,y,1,1))
163     !
164     status = nf90_inq_varid(ncid, "TRNDIC", VarId) !PISCES
165     IF (status == nf90_noerr) THEN
166        CALL Read_Ncdf_var('TRNDIC',TRIM(restart_trc_file),tabvar1,1,k)
167     ELSE
168        status = nf90_inq_varid(ncid, "TRNNO3"  , VarId) ! LOBSTER
169        IF (status == nf90_noerr) THEN
170           CALL Read_Ncdf_var('TRNNO3',TRIM(restart_trc_file),tabvar1,1,k)
171        ELSE
172           status = nf90_inq_varid(ncid, "TRNCFC11", VarId) ! CFC
173           IF (status == nf90_noerr) THEN
174              CALL Read_Ncdf_var('TRNCFC11',TRIM(restart_trc_file),tabvar1,1,k)
175           ELSE
176              status = nf90_inq_varid(ncid, "TRNCLR  ", VarId) ! My TRC
177              IF (status == nf90_noerr) THEN
178                 CALL Read_Ncdf_var('TRNCLR',TRIM(restart_trc_file),tabvar1,1,k)
179              ELSE
180                 WRITE(*,*) 'No suitable tracer found to build the mask '
181              ENDIF
182           ENDIF
183        ENDIF
184     ENDIF
185     WHERE( tabvar1(:,:,1,1) == 0. ) 
186        G0%tmask(:,:,k) = 0.
187     END WHERE
188     DEALLOCATE(tabvar1)
189  END DO
190
191  !
192  ! write dimensions in output file
193  !         
194  CALL Write_Ncdf_dim('x',Child_file,nxfin)
195  CALL Write_Ncdf_dim('y',Child_file,nyfin)
196  CALL Write_Ncdf_dim('z',Child_file,z)
197  CALL Write_Ncdf_dim(TRIM(timedimname),Child_file,0) 
198  !
199  !
200  VARIABLE :  DO i = 1,SIZE(Ncdf_varname)     
201     !     
202     ! loop on variables names
203     !     
204     SELECT CASE (TRIM(Ncdf_varname(i)))
205        !
206        !copy nav_lon from child coordinates to output file     
207        !
208     CASE('nav_lon')
209        WRITE(*,*) 'copy nav_lon'
210        CALL Read_Ncdf_var('nav_lon',TRIM(Childcoordinates),tabtemp2D) 
211        CALL Write_Ncdf_var('nav_lon',(/'x','y'/),Child_file,tabtemp2D,'float')
212        CALL Copy_Ncdf_att('nav_lon',TRIM(restart_trc_file),Child_file, &
213             MINVAL(tabtemp2D),MAXVAL(tabtemp2D))
214        DEALLOCATE(tabtemp2D)
215        Interpolation = .FALSE.
216        !       
217        !copy nav_lat from child coordinates to output file
218        !
219     CASE('nav_lat')             
220        WRITE(*,*) 'copy nav_lat'
221        CALL Read_Ncdf_var('nav_lat',TRIM(Childcoordinates),tabtemp2D) 
222        CALL Write_Ncdf_var('nav_lat',(/'x','y'/),Child_file,tabtemp2D,'float')
223        CALL Copy_Ncdf_att('nav_lat',TRIM(restart_trc_file),Child_file, &
224             MINVAL(tabtemp2D),MAXVAL(tabtemp2D)) 
225        DEALLOCATE(tabtemp2D)
226        Interpolation = .FALSE.
227        !
228        !copy nav_lev from restart_file to output file
229        !
230     CASE('nav_lev')
231
232        WRITE(*,*) 'copy nav_lev'
233        CALL Read_Ncdf_var('nav_lev',TRIM(restart_trc_file),nav_lev) 
234        CALL Write_Ncdf_var('nav_lev','z',Child_file,nav_lev,'float')
235        CALL Copy_Ncdf_att('nav_lev',TRIM(restart_trc_file),Child_file)     
236        Interpolation = .FALSE.
237        !
238        !copy time from restart_file to output file                       
239        !
240     CASE('time_counter')
241        WRITE(*,*) 'copy time_counter'
242        CALL Read_Ncdf_var('time_counter',TRIM(restart_trc_file),tabtemp1D) 
243        tabtemp1D = tabtemp1D * rhot
244        CALL Write_Ncdf_var('time_counter',TRIM(timedimname),Child_file,tabtemp1D,'double')
245        CALL Copy_Ncdf_att('time_counter',TRIM(restart_trc_file),Child_file) 
246        DEALLOCATE(tabtemp1D)
247        Interpolation = .FALSE.
248        !
249        !copy info from restart_file to output file
250        !
251     CASE('arak0') 
252        WRITE(*,*) 'copy trp info'       
253        CALL Read_Ncdf_var('arak0',TRIM(restart_trc_file),tabtemp0dreal) 
254        CALL Write_Ncdf_var('arak0',Child_file,tabtemp0dreal,'double')
255        CALL Copy_Ncdf_att('arak0',TRIM(restart_trc_file),Child_file) 
256        Interpolation = .FALSE. 
257        !
258     CASE('kt') 
259        WRITE(*,*) 'copy kt'
260        CALL Read_Ncdf_var('kt',TRIM(restart_trc_file),tabtemp0dreal) 
261        tabtemp0dreal = tabtemp0dreal * rhot
262        CALL Write_Ncdf_var('kt',Child_file,tabtemp0dreal,'double')
263        CALL Copy_Ncdf_att('kt',TRIM(restart_trc_file),Child_file) 
264        Interpolation = .FALSE.
265        !
266     CASE DEFAULT
267        varname = Ncdf_varname(i) 
268        WRITE(*,*) TRIM(varname),' interpolation ...'     
269        vert_coord_name = 'z'
270        posvar='T'
271        Interpolation = .TRUE.           
272        !     
273     END SELECT
274     IF( Interpolation ) THEN
275        !       
276        nbvert_lev = z 
277        !                 
278        t = 1
279
280        ALLOCATE(detected_pts(SIZE(G0%tmask,1),SIZE(G0%tmask,2),nbvert_lev))                             
281        ALLOCATE(tabvar1(x,y,1,2))
282        ALLOCATE(tabvar2(x,y,1,1))
283        ALLOCATE(tabvar3(x,y,1,1))
284        ALLOCATE(masksrc(x,y))
285        ALLOCATE(tabinterp4d(nxfin,nyfin,1,1)) 
286
287        !       
288        DO n = 1,nbvert_lev
289           !
290           WRITE(*,*) 'interpolate/extrapolate ', &
291                TRIM(varname),' for vertical level = ',n   
292           !
293           CALL Read_Ncdf_var(varname,TRIM(restart_trc_file),tabvar1,t,n)
294           IF(n==1) THEN
295              !                           
296           ELSE IF (n==2) THEN
297              tabvar2(:,:,:,1) = tabvar1(:,:,:,2)
298           ELSE
299              tabvar3(:,:,:,1) = tabvar2(:,:,:,1)
300              tabvar2(:,:,:,1) = tabvar1(:,:,:,2)
301           ENDIF
302           !
303           !
304           IF(MAXVAL(G1%tmask(:,:,n)) == 0.) THEN
305              tabinterp4d = 0.0
306              WRITE(*,*) 'only land points on level ',n                     
307           ELSE
308              CALL extrap_detect(G0,G1,detected_pts(:,:,n),n)                                           
309
310              CALL correct_field(detected_pts(:,:,n),tabvar1,tabvar2,&
311                   tabvar3,G0,nav_lev,masksrc,n)                                 
312
313              SELECT CASE(TRIM(interp_type))
314              CASE('bilinear')                                                       
315                 CALL get_remap_matrix(G0%nav_lat,G1%nav_lat, &
316                      G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
317                 CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
318                      matrix,src_add,dst_add)     
319              CASE('bicubic')                                   
320                 CALL get_remap_bicub(G0%nav_lat,G1%nav_lat, &
321                      G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
322                 CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),&
323                      nxfin,nyfin,matrix,src_add,dst_add) 
324              END SELECT
325              !                     
326              CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
327                   G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin-jpizoom+1,jmin-jpjzoom+1)
328           ENDIF
329
330           tabinterp4d(:,:,1,1) =  tabinterp4d(:,:,1,1) * G1%tmask(:,:,n)
331           ztrmas = 0.0
332           ztrcor = 0.0
333           DO jj=1, nyfin
334           DO ji=1, nxfin
335              zvolk = cvol(ji,jj,n)
336              ztrcor = ztrcor + MIN( 0., tabinterp4d(ji,jj,1,1) ) * zvolk
337              tabinterp4d(ji,jj,1,1) = MAX( 0., tabinterp4d(ji,jj,1,1)) 
338              ztrmas  = ztrmas + tabinterp4d(ji,jj,1,1)  * zvolk
339           END DO
340           END DO
341           IF( ztrcor .LT. 0. ) THEN
342               WRITE(*,*) 'Correcting negative concentration'   
343               zcoef = 1. + ztrcor / ztrmas
344               tabinterp4d(:,:,1,1) =  tabinterp4d(:,:,1,1) *zcoef * G1%tmask(:,:,n)               
345           ENDIF       
346           IF(MINVAL(tabinterp4d(:,:,1,1)) .LT.0. ) STOP
347
348           !
349
350           status = nf90_inq_varid(ncid,TRIM(varname) , VarId)
351           status = nf90_inquire_variable(ncid, VarId, ndims=numdims)
352           IF( numdims == 3) THEN
353              dimnames(1)='x'
354              dimnames(2)='y'
355              dimnames(3)=TRIM(timedimname)
356              ALLOCATE(tabvar3d(SIZE(tabinterp4d,1),SIZE(tabinterp4d,2),SIZE(tabinterp4d,3)))
357              tabvar3d=tabinterp4d(:,:,:,1)
358              CALL Write_Ncdf_var(TRIM(varname),dimnames, &
359                   Child_file,tabvar3d,t,'double')
360              DEALLOCATE(tabvar3d)
361           ELSE       
362              dimnames(1)='x'
363              dimnames(2)='y'
364              dimnames(3)=vert_coord_name
365              dimnames(4)=TRIM(timedimname)
366              CALL Write_Ncdf_var(TRIM(varname),dimnames, &
367                   Child_file,tabinterp4d,t,n,'double')
368           ENDIF
369           !
370           CALL Copy_Ncdf_att(TRIM(varname),TRIM(restart_trc_file),Child_file)
371           !
372           IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,src_add,dst_add)   
373           !
374           IF( numdims == 3) CYCLE VARIABLE
375           !
376        END DO
377        !
378        DEALLOCATE(detected_pts)
379        DEALLOCATE(tabinterp4d)
380        DEALLOCATE(tabvar1,tabvar2,tabvar3)
381        DEALLOCATE(masksrc) 
382        !                     
383     ENDIF
384
385     prefix = varname(1:3)
386     suffix = varname(4:LEN_TRIM(varname))   
387
388     IF(rhot == 1 .OR. prefix/= 'TRB') THEN
389        WRITE(*,*) ''   
390        WRITE(*,*) 'no time interpolation for', varname
391     ELSE   
392        WRITE(*,*) ''
393        WRITE(*,*) 'time interpolation for', varname
394        ALLOCATE(trn(nxfin,nyfin,z,1),trb(nxfin,nyfin,z,1))
395        varname2 = 'TRN'//TRIM(suffix)
396        now_wght = (rhot-1.)/rhot
397        before_wght = 1./rhot
398        !   
399        CALL Read_Ncdf_var(TRIM(varname),Child_file,trb)
400        CALL Read_Ncdf_var(TRIM(varname2),Child_file,trn)           
401        trb = now_wght*trn + before_wght*trb
402
403        dimnames(1)='x'
404        dimnames(2)='y'
405        dimnames(3)='z'
406        dimnames(4)=TRIM(timedimname)
407        CALL Write_Ncdf_var(TRIM(varname),dimnames,Child_file,trb,'double')
408        !   
409     ENDIF
410
411  END DO VARIABLE
412  !
413  WRITE(*,*) ' '
414  WRITE(*,*) '******* restart file successfully created *******' 
415  WRITE(*,*) ' '
416  !
417  STOP
418END PROGRAM create_rstrt_trc
419
420
Note: See TracBrowser for help on using the repository browser.