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_trc.f90 in utils/tools/NESTING/src – NEMO

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

Last change on this file since 12253 was 10398, checked in by clem, 5 years 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
RevLine 
[2136]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       
[10398]19  USE agrif_connect_topo
[2136]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 
[10398]43  CHARACTER*100 :: Child_file,Childcoordinates,varname,varname2,Child_Bathy_Level,Child_Bathy_Meter   
44  REAL*8, POINTER, DIMENSION(:,:,:) :: tabvar3d => NULL()
[2136]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()
[10398]48  REAL*8, POINTER, DIMENSION(:,:) :: tabtemp2D => NULL()
[2136]49  INTEGER,DIMENSION(:),POINTER :: src_add,dst_add  => NULL()
[10398]50  REAL*8,DIMENSION(:,:),POINTER :: matrix => NULL()
[2136]51  LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL()
52  LOGICAL, DIMENSION(:,:,:), POINTER :: detected_pts
[10398]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
[2136]57  CHARACTER(len=20),DIMENSION(4) :: dimnames
58  CHARACTER(len=80) :: namelistname
59  TYPE(Coordinates) :: G0,G1
60  REAL*8 :: tabtemp0dreal
61
[10398]62  LOGICAL, PARAMETER :: conservation = .FALSE.
[2136]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  !     
[10381]75  IF(TRIM(restart_trc_file) == '') THEN
[2136]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)   
[10398]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)   
[2136]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) 
[10398]103  CALL Read_Ncdf_dim('nav_lev',restart_trc_file,z)
[2136]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  !
[10398]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
[10381]141  CALL Init_tmask(parent_bathy_level,G0,x,y)
[10398]142  CALL Init_tmask(child_bathy_level,G1,nxfin,nyfin)
[2136]143 
[10398]144!!clem  G0%tmask = 1.   
[2136]145
[10398]146  ! which dataset
[2136]147  status = nf90_open(TRIM(restart_trc_file), NF90_NOWRITE, ncid) ! Open dataset
[10398]148  DO jk = 1, z
[2136]149     ALLOCATE(tabvar1(x,y,1,1))
150     !
151     status = nf90_inq_varid(ncid, "TRNDIC", VarId) !PISCES
152     IF (status == nf90_noerr) THEN
[10398]153        CALL Read_Ncdf_var('TRNDIC',TRIM(restart_trc_file),tabvar1,1,jk)
[2136]154     ELSE
155        status = nf90_inq_varid(ncid, "TRNNO3"  , VarId) ! LOBSTER
156        IF (status == nf90_noerr) THEN
[10398]157           CALL Read_Ncdf_var('TRNNO3',TRIM(restart_trc_file),tabvar1,1,jk)
[2136]158        ELSE
159           status = nf90_inq_varid(ncid, "TRNCFC11", VarId) ! CFC
160           IF (status == nf90_noerr) THEN
[10398]161              CALL Read_Ncdf_var('TRNCFC11',TRIM(restart_trc_file),tabvar1,1,jk)
[2136]162           ELSE
163              status = nf90_inq_varid(ncid, "TRNCLR  ", VarId) ! My TRC
164              IF (status == nf90_noerr) THEN
[10398]165                 CALL Read_Ncdf_var('TRNCLR',TRIM(restart_trc_file),tabvar1,1,jk)
[2136]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. ) 
[10398]173        G0%tmask(:,:,jk) = 0.
[2136]174     END WHERE
175     DEALLOCATE(tabvar1)
176  END DO
177
178  !
179  ! write dimensions in output file
[10398]180  WRITE(*,*) 'write dimensions'
[2136]181  !         
182  CALL Write_Ncdf_dim('x',Child_file,nxfin)
183  CALL Write_Ncdf_dim('y',Child_file,nyfin)
[10398]184  CALL Write_Ncdf_dim('nav_lev',Child_file,z)
185  CALL Write_Ncdf_dim('time_counter',Child_file,0) 
[2136]186  !
187  !
[10398]188  DO ii = 1,SIZE(Ncdf_varname)     
[2136]189     !     
190     ! loop on variables names
[10398]191     varname = TRIM(Ncdf_varname(ii))
192     WRITE(*,*) 'var = ',TRIM(varname)     
[2136]193     !     
[10398]194     SELECT CASE (TRIM(varname))
[2136]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')
[10398]199        CALL Copy_Ncdf_att('nav_lon',TRIM(restart_trc_file),Child_file,MINVAL(tabtemp2D),MAXVAL(tabtemp2D))
[2136]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')
[10398]206        CALL Copy_Ncdf_att('nav_lat',TRIM(restart_trc_file),Child_file,MINVAL(tabtemp2D),MAXVAL(tabtemp2D)) 
[2136]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
[10398]219        CALL Write_Ncdf_var('time_counter','time_counter',Child_file,tabtemp1D,'double')
[2136]220        CALL Copy_Ncdf_att('time_counter',TRIM(restart_trc_file),Child_file) 
221        DEALLOCATE(tabtemp1D)
222        Interpolation = .FALSE.
223        !
[10398]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) 
[2136]234        Interpolation = .FALSE.
235        !
236     CASE DEFAULT
[10398]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
[2136]242        posvar='T'
243        Interpolation = .TRUE.           
244        !     
245     END SELECT
[10398]246     !     
247     ! --- start interpolation --- !
[2136]248     IF( Interpolation ) THEN
249        !       
[10398]250        IF( vert_coord_name == '1' ) THEN
251           nbvert_lev = 1
252        ELSE
253           nbvert_lev = z
254        ENDIF
[2136]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           !
[10398]266           WRITE(*,*) 'interpolate/extrapolate for vertical level = ',n   
[2136]267           !
[10398]268           CALL Read_Ncdf_var(varname,TRIM(restart_trc_file),tabvar1,1,n)
[2136]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
[10398]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             
[2136]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              !                     
[10398]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
[2136]309           ENDIF
310
[10398]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
[2136]316
317           !
[10398]318           dimnames(1)='x'
319           dimnames(2)='y'
320           IF( vert_coord_name == '1' ) THEN
321              dimnames(3)='time_counter'
322             
[2136]323              ALLOCATE(tabvar3d(SIZE(tabinterp4d,1),SIZE(tabinterp4d,2),SIZE(tabinterp4d,3)))
324              tabvar3d=tabinterp4d(:,:,:,1)
[10398]325              CALL Write_Ncdf_var(TRIM(varname),dimnames,Child_file,tabvar3d,1,'double')
[2136]326              DEALLOCATE(tabvar3d)
[10398]327           ELSE
[2136]328              dimnames(3)=vert_coord_name
[10398]329              dimnames(4)='time_counter'
330             
331              CALL Write_Ncdf_var(TRIM(varname),dimnames,Child_file,tabinterp4d,1,n,'double')
[2136]332           ENDIF
[10398]333           !             
[2136]334           !
335           CALL Copy_Ncdf_att(TRIM(varname),TRIM(restart_trc_file),Child_file)
336           !
[10398]337           !
[2136]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
[10398]350     ! change the before fields
[2136]351     prefix = varname(1:3)
352     suffix = varname(4:LEN_TRIM(varname))   
353
354     IF(rhot == 1 .OR. prefix/= 'TRB') THEN
355        WRITE(*,*) ''   
[10398]356        WRITE(*,*) 'no time interpolation for ',TRIM(varname)
[2136]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'
[10398]369        dimnames(3)='nav_lev'
370        dimnames(4)='time_counter'
[2136]371        CALL Write_Ncdf_var(TRIM(varname),dimnames,Child_file,trb,'double')
[10398]372        DEALLOCATE(trn,trb)
[2136]373        !   
374     ENDIF
375
[10398]376  END DO
[2136]377  !
378  WRITE(*,*) ' '
[10398]379  WRITE(*,*) ' --- list of all the variables that have been interpolated --- '
380  WRITE(*,*) Ncdf_varname
381  WRITE(*,*) ' '
[2136]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.