Ignore:
Timestamp:
2018-12-11T21:00:08+01:00 (2 years ago)
Author:
clem
Message:

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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • utils/tools/NESTING/src/agrif_create_bathy.f90

    r10093 r10381  
    4141  !       
    4242  CHARACTER(len=80) ::   namelistname 
    43   CHARACTER*100     ::   Childmeter_file,Childlevel_file,Child_coordinates,child_ps, Child_domcfg     
     43  CHARACTER*100     ::   child_coordinates, child_level, child_meter, child_domcfg     
    4444  LOGICAL ::   identical_grids      
    4545  INTEGER ::   nbadd,status,narg,iargc      
     
    6969  ! read input file (namelist.input) 
    7070  CALL read_namelist(namelistname) 
    71   !       
     71 
     72  ! if level or meter name is missing 
     73  IF( TRIM(parent_level_name) == '' )   parent_level_name='mbathy' 
     74  IF( TRIM(parent_meter_name) == '' )   parent_meter_name='Bathymetry' 
     75 
    7276  ! define names of child grid files 
    73   CALL set_child_name(parent_coordinate_file,Child_coordinates)  
    74   IF( TRIM(parent_bathy_file) .NE. '/NULL' )    CALL set_child_name(parent_bathy_file,Childlevel_file)             
    75   IF( TRIM(parent_bathy_file) .NE. '/NULL' )    CALL set_child_name(parent_domcfg_output,Child_domcfg)  
     77  CALL set_child_name(parent_coordinate_file,child_coordinates)  
     78  IF( TRIM(parent_bathy_level) /= '' )   CALL set_child_name(parent_bathy_level,child_level)             
     79  IF( TRIM(parent_bathy_meter) /= '' )   CALL set_child_name(parent_bathy_meter,child_meter) 
     80  IF( TRIM(parent_domcfg_out)  /= '' )   CALL set_child_name(parent_domcfg_out,child_domcfg)  
    7681  ! 
     82  IF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) == '') THEN 
     83     WRITE(*,*) 'ERROR ***** one needs at least to define parent_bathy_level or parent_bathy_meter ...' 
     84     STOP 
     85  ENDIF 
    7786  !                                                   ! ------------------------------------------------------------------ 
    7887  IF( .NOT.new_topo .AND. .NOT.partial_steps ) THEN   ! **** First option : no new topo file & no partial steps 
     
    8190     WRITE(*,*) '*** First option: no new topo file & no partial steps' 
    8291     !       
    83      ! read coarse grid bathymetry and coordinates 
     92     ! read coarse grid coordinates 
    8493     status = Read_Coordinates(TRIM(parent_coordinate_file),G0)     
    85      status = Read_bathy_level(TRIM(parent_bathy_file),G0) 
     94 
     95     ! read coarse grid bathymetry 
     96     IF( TRIM(parent_bathy_level) /= '' ) THEN 
     97        status = Read_bathy_level(TRIM(parent_bathy_level),G0) 
     98        CALL levels_to_meter(G0) 
     99     ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN 
     100        status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 
     101        CALL meter_to_levels(G0) 
     102     ENDIF 
    86103     !            
    87104     ! read fine grid coordinates 
    88      status = Read_Coordinates(TRIM(Child_coordinates),G1,pacifique) 
     105     status = Read_Coordinates(TRIM(child_coordinates),G1,pacifique) 
    89106     ! 
    90107     ! stop if error in size 
     
    113130        WHERE(G0%nav_lon < 0.001)   G0%nav_lon = G0%nav_lon + 360. 
    114131     ENDIF 
    115  
    116      ! From levels to meters 
     132     !              
    117133     ALLOCATE(G1%bathy_meter(nxfin,nyfin)) 
    118      CALL levels_to_meter(G0) 
    119      !              
     134     ! 
    120135     ! compute remapping matrix thanks to SCRIP package 
    121136     CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add) 
     
    189204        WHERE (bathy_test.EQ.0.)   G1%bathy_level = 0. 
    190205 
    191         DEALLOCATE(bathy_test)            
    192      ENDIF 
     206        DEALLOCATE(bathy_test) 
     207     ENDIF 
     208     ! 
     209     CALL levels_to_meter(G1) ! needed for domcfg 
    193210     ! 
    194211     ! store interpolation result in output file 
    195      CALL levels_to_meter(G1)   ! From levels to meters 
    196      status = Write_Bathy_level(TRIM(Childlevel_file),G1) 
    197      status = write_domcfg(TRIM(Child_domcfg),G1) 
     212     status = Write_Bathy_level(TRIM(child_level),G1) 
     213     status = write_domcfg(TRIM(child_domcfg),G1) 
     214     !!status = write_domcfg(TRIM(parent_domcfg_out),G0) ! do not activate it 
    198215 
    199216     WRITE(*,*) '****** Bathymetry successfully created if no new topo ******' 
     
    203220  ELSE                                                ! **** Second option : new topo file or partial steps      
    204221     !                                                ! -----------------------------------------------------  
    205      WRITE(*,*) '***Second option : new topo or partial steps' 
     222     WRITE(*,*) '*** Second option : new topo or partial steps' 
    206223 
    207224     ! read fine and coarse grids coordinates file 
    208225     status = Read_Coordinates(TRIM(parent_coordinate_file),G0) 
    209      status = Read_Coordinates(TRIM(Child_coordinates),G1,Pacifique) 
     226     status = Read_Coordinates(TRIM(child_coordinates),G1,Pacifique) 
    210227     ! 
    211228     ! check error in size 
     
    221238     ENDIF 
    222239     ! 
    223      ! === From here on: G0 is the grid associated with the new topography (like gebco or etopo) === 
     240     ! === From here on: G0 is the grid associated with the new topography (as gebco or etopo) === 
    224241     ! 
    225242     ! read bathymetry data set => G0%bathy_meter 
    226243     IF( new_topo ) THEN                                       ! read G0%bathy_meter (on a reduced grid) and G1 coordinates 
    227244        DEALLOCATE( G0%nav_lon, G0%nav_lat ) 
    228         status = Read_bathy_meter(TRIM(elevation_database),G0,G1,Pacifique) 
     245        status = read_bathy_coord(TRIM(elevation_database),G0,G1,Pacifique) 
    229246     ELSE                                                      ! read G0%bathy_meter (on the global grid) 
    230         status = Read_Bathymeter(TRIM(parent_bathy_meter),G0)    
     247        IF( TRIM(parent_bathy_level) /= '' ) THEN 
     248           status = Read_bathy_level(TRIM(parent_bathy_level),G0) 
     249           CALL levels_to_meter(G0) 
     250        ELSEIF( TRIM(parent_bathy_level) == '' .AND. TRIM(parent_bathy_meter) /= '') THEN 
     251           status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 
     252        ENDIF 
    231253        IF(Pacifique) THEN 
    232254           WHERE(G0%nav_lon < 0.0001)   G0%nav_lon = G0%nav_lon + 360. 
     
    408430     ! 
    409431     ! Coarse grid bathymetry : G0%bathy_meter (on the global grid) 
    410      IF( partial_steps ) THEN                      
    411         status = Read_Bathymeter(TRIM(parent_bathy_meter),G0) 
     432     IF( TRIM(parent_bathy_meter) /= '') THEN 
     433        status = read_bathy_meter(TRIM(parent_bathy_meter),G0) 
    412434     ELSE 
    413         status = Read_Bathymeter(TRIM(parent_bathy_file),G0) 
     435        status = Read_bathy_level(TRIM(parent_bathy_level),G0) 
     436        CALL levels_to_meter(G0) 
    414437     ENDIF 
    415438      
     
    586609     ENDIF 
    587610     ! 
    588      !                              ! ---------------- 
    589      IF( partial_steps ) THEN       ! If partial steps 
    590         !                           ! ---------------- 
    591         ! 
    592         CALL get_partial_steps(G1)  ! recompute gdepw_ps for G1 
    593  
    594         IF(bathy_update)   CALL Update_Parent_Bathy( G0,G1 )  
    595         CALL set_child_name(parent_bathy_meter,child_ps) 
    596          
    597                            status = Write_Bathy_meter(TRIM(child_ps),G1) 
    598         IF(bathy_update)   status = Write_Bathy_meter(TRIM(updated_parent_file),G0) 
    599         IF(bathy_update)   status = write_domcfg(TRIM(updated_parent_domcfg),G0) 
    600  
    601         CALL get_partial_steps(G1) 
    602         ! 
    603         G1%bathy_level=NINT(G1%bathy_level) 
    604         ! 
    605         ! store interpolation result in output file 
    606         CALL levels_to_meter(G1)   ! From levels to meters 
    607         IF( TRIM(parent_bathy_file) .NE. '/NULL' )    status = Write_Bathy_level(TRIM(Childlevel_file),G1) 
    608         IF( TRIM(parent_bathy_file) .NE. '/NULL' )    status = write_domcfg(TRIM(Child_domcfg),G1) 
    609         ! 
    610         WRITE(*,*) '****** Bathymetry successfully created for partial cells ******' 
    611         STOP 
    612         ! 
    613         !                           ! ---------------- 
    614      ELSE                           ! No partial steps 
    615         !                           ! ---------------- 
    616         ! 
    617         ! convert bathymetry from meters to levels 
    618         CALL meter_to_levels(G1)  
    619  
    620         IF(bathy_update)   CALL Update_Parent_Bathy( G0,G1 )  
    621         CALL set_child_name(parent_bathy_meter,child_ps) 
    622          
    623                            status = Write_Bathy_meter(TRIM(child_ps),G1)         
    624         IF(bathy_update)   status = Write_Bathy_meter(TRIM(updated_parent_file),G0) 
    625         IF(bathy_update)   status = write_domcfg(TRIM(updated_parent_domcfg),G0) 
    626         ! 
    627         ! store interpolation result in output file 
    628         CALL levels_to_meter(G1)   ! From levels to meters 
    629         status = Write_Bathy_level(TRIM(Childlevel_file),G1) 
    630         status = write_domcfg(TRIM(Child_domcfg),G1) 
    631         ! 
    632         WRITE(*,*) '****** Bathymetry successfully created for full cells ******' 
    633         STOP 
    634         !  
    635      ENDIF 
     611     IF( partial_steps ) THEN 
     612        CALL get_partial_steps(G1)  ! recompute bathy_level and gdepw_ps for G1 (and correct bathy_meter) 
     613     ELSE 
     614        CALL meter_to_levels(G1)    ! convert bathymetry from meters to levels 
     615     ENDIF 
     616     ! 
     617     ! update parent grid 
     618     IF(bathy_update)   CALL Update_Parent_Bathy( G0,G1 )  
     619      
     620     IF(bathy_update)   status = Write_Bathy_meter(TRIM(parent_bathy_meter_updated),G0) 
     621     IF(bathy_update)   status = write_domcfg(TRIM(parent_domcfg_updated),G0) 
     622      
     623     ! store interpolation result in output file 
     624     IF( TRIM(parent_bathy_level) /= '' )   status = Write_Bathy_level(TRIM(child_level),G1) 
     625     IF( TRIM(parent_bathy_meter) /= '' )   status = Write_Bathy_meter(TRIM(child_meter),G1) 
     626     IF( TRIM(parent_domcfg_out)  /= '' )   status = write_domcfg(TRIM(child_domcfg),G1) 
     627     ! 
     628     WRITE(*,*) '****** Bathymetry successfully created ******' 
     629     STOP 
    636630  ENDIF 
    637631  ! 
Note: See TracChangeset for help on using the changeset viewer.