source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/TOOLS/NESTING/src/agrif_create_restart_trc.f90 @ 5445

Last change on this file since 5445 was 5445, checked in by davestorkey, 5 years ago

Clear SVN keywords from 2015/dev_r5021_UKMO1_CICE_coupling branch.

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