Changeset 882


Ignore:
Timestamp:
06/11/19 14:59:17 (5 years ago)
Author:
ymipsl
Message:

Metric is now write in start.nc/restart.nc
Metric can be read at restart if read_metric=y.

YM

Location:
codes/icosagcm/trunk/src
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/icosa_init.f90

    r710 r882  
    2222  USE diagflux_mod 
    2323  USE profiling_mod 
     24  USE read_metric_mod 
     25  USE xios 
    2426  IMPLICIT NONE 
    2527   
    … …  
    3941  !$OMP PARALLEL   
    4042    CALL switch_omp_no_distrib_level 
     43    CALL read_metric 
    4144    CALL compute_geometry 
    4245    CALL check_total_area 
    … …  
    6467    CALL init_diagflux 
    6568    CALL zero_du_phys 
     69     
     70    IF (is_omp_master) CALL xios_timer_resume("dynamico") 
    6671    CALL timeloop 
     72    IF (is_omp_master) CALL xios_timer_suspend("dynamico") 
     73     
    6774    CALL switch_omp_no_distrib_level 
    6875  !$OMP END PARALLEL 

  • codes/icosagcm/trunk/src/output/write_etat0.f90

    r581 r882  
    2626    TYPE(t_field),POINTER,SAVE :: f_ulat(:) 
    2727    TYPE(t_field),POINTER,SAVE :: f_theta_rhodz_1d(:) 
     28    TYPE(t_field),POINTER,SAVE :: f_xcell(:),f_ycell(:),f_zcell(:) 
    2829    REAL(rstd), POINTER :: theta_rhodz(:,:,:),theta_rhodz_1d(:,:) 
     30    REAL(rstd), POINTER :: xcell(:), ycell(:), zcell(:) 
    2931    INTEGER :: ind 
    3032     
    … …  
    3335    CALL allocate_field(f_ulat,field_t,type_real,llm,name='ulat') 
    3436    CALL allocate_field(f_theta_rhodz_1d,field_t,type_real,llm,name='theta_rhodz') 
     37    CALL allocate_field(f_xcell,field_t,type_real,name='xcell') 
     38    CALL allocate_field(f_ycell,field_t,type_real,name='ycell') 
     39    CALL allocate_field(f_zcell,field_t,type_real,name='zcell') 
    3540 
    3641!$OMP BARRIER     
    3742    DO ind=1, ndomain 
    3843       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
     44       CALL swap_dimensions(ind) 
     45       CALL swap_geometry(ind) 
    3946       theta_rhodz=f_theta_rhodz(ind) ; theta_rhodz_1d=f_theta_rhodz_1d(ind) 
    4047       theta_rhodz_1d(:,:)=theta_rhodz(:,:,1) 
     48       xcell=f_xcell(ind) ; xcell=xyz_i(:,1)/radius 
     49       ycell=f_ycell(ind) ; ycell=xyz_i(:,2)/radius 
     50       zcell=f_zcell(ind) ; zcell=xyz_i(:,3)/radius 
    4151    ENDDO 
    4252     
    … …  
    4555 
    4656    IF(hydrostatic) THEN 
    47        CALL write_restart(it,f_ps,f_phis,f_theta_rhodz_1d,f_u, f_ulon, f_ulat, f_q) 
     57       CALL write_restart(it,f_ps,f_phis,f_theta_rhodz_1d,f_u, f_ulon, f_ulat, f_q, f_xcell, f_ycell, f_zcell ) 
    4858    ELSE 
    49        CALL write_restart(it,f_ps,f_phis,f_theta_rhodz_1d,f_u, f_ulon, f_ulat, f_q, f_geopot, f_W) 
     59       CALL write_restart(it,f_ps,f_phis,f_theta_rhodz_1d,f_u, f_ulon, f_ulat, f_q, f_geopot, f_W, f_xcell, f_ycell, f_zcell) 
    5060    END IF 
    5161    CALL deallocate_field(f_ulon) 

  • codes/icosagcm/trunk/src/output/xios_mod.F90

    r707 r882  
    3131     using_xios=.TRUE. 
    3232     CALL xios_context_initialize("icosagcm",comm_icosa) 
    33      CALL xios_get_handle("icosagcm",ctx_hdl) 
    34      CALL xios_set_current_context(ctx_hdl) 
     33     CALL xios_context_initialize("icosagcm_input",comm_icosa) 
     34     CALL xios_set_context 
    3535    
    3636 END SUBROUTINE xios_init   
    … …  
    6060  TYPE(t_domain),POINTER :: d 
    6161 
     62   CALL xios_set_context 
    6263!$OMP BARRIER 
    6364!$OMP MASTER 
    6465!   CALL xios_context_initialize("icosagcm",comm_icosa) 
    65    CALL xios_get_handle("icosagcm",ctx_hdl) 
    66    CALL xios_set_current_context(ctx_hdl) 
     66!   CALL xios_get_handle("icosagcm",ctx_hdl) 
     67!   CALL xios_set_current_context(ctx_hdl) 
    6768   lev_value(:) = (/ (l,l=1,llm) /)      
    6869   lev_valuep1(:) = (/ (l,l=1,llm+1) /)      
    … …  
    274275    
    275276 END SUBROUTINE xios_init_write_field 
    276    
     277  
     278  
     279 SUBROUTINE xios_init_write_field_input 
     280 USE genmod 
     281 USE mpipara 
     282 USE xios 
     283 USE grid_param 
     284 USE domain_mod 
     285 USE dimensions 
     286 USE spherical_geom_mod 
     287 USE geometry 
     288 USE mpi_mod 
     289 USE time_mod 
     290 USE metric, ONLY : vup,vdown, cell_glo 
     291 USE icosa,ONLY  : getin 
     292 IMPLICIT NONE 
     293  TYPE(xios_context) :: ctx_hdl 
     294  TYPE(xios_duration)      :: dtime 
     295  INTEGER :: ncell, ncell_tot, ncell_glo(0:mpi_size-1), displ 
     296  INTEGER :: ind, i,j,k,l,ij 
     297  REAL(rstd),ALLOCATABLE    :: lon(:), lat(:), bounds_lon(:,:), bounds_lat(:,:) 
     298  INTEGER, ALLOCATABLE      :: ind_glo(:) 
     299  TYPE(t_domain),POINTER :: d 
     300  CHARACTER(len=255) :: etat0_type 
     301  LOGICAL :: read_metric_ 
     302 
     303   CALL xios_set_context_input 
     304!$OMP BARRIER 
     305!$OMP MASTER 
     306    
     307   ncell=0 
     308   DO ind=1,ndomain 
     309     d=>domain(ind) 
     310         
     311     DO j=d%jj_begin,d%jj_end 
     312       DO i=d%ii_begin,d%ii_end 
     313         IF (domain(ind)%own(i,j)) ncell=ncell+1 
     314       ENDDO 
     315     ENDDO 
     316   ENDDO       
     317   ncell_i=ncell 
     318    
     319   CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) 
     320 
     321   displ=0 
     322   DO i=1,mpi_rank 
     323     displ=displ+ncell_glo(i-1) 
     324   ENDDO 
     325 
     326   ncell_tot=sum(ncell_glo(:)) 
     327    
     328   ALLOCATE(lon(ncell), lat(ncell), bounds_lon(0:5,ncell), bounds_lat(0:5,ncell), ind_glo(ncell))  
     329    
     330   ncell=0 
     331   DO ind=1,ndomain 
     332     d=>domain(ind) 
     333         
     334     DO j=d%jj_begin,d%jj_end 
     335       DO i=d%ii_begin,d%ii_end 
     336         IF (domain(ind)%own(i,j)) THEN  
     337           ncell=ncell+1 
     338           CALL xyz2lonlat(d%xyz(:,i,j),lon(ncell),lat(ncell)) 
     339           lon(ncell)=lon(ncell)*180/Pi 
     340           lat(ncell)=lat(ncell)*180/Pi 
     341           DO k=0,5 
     342             CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,ncell), bounds_lat(k,ncell)) 
     343             bounds_lat(k,ncell)=bounds_lat(k,ncell)*180/Pi 
     344             bounds_lon(k,ncell)=bounds_lon(k,ncell)*180/Pi 
     345           ENDDO 
     346           ind_glo(ncell)=domain(ind)%assign_cell_glo(i,j)-1  
     347         ENDIF 
     348       ENDDO 
     349     ENDDO 
     350   ENDDO          
     351 
     352   CALL xios_set_domain_attr("i",ni_glo=ncell_tot, ibegin=displ, ni=ncell) 
     353   CALL xios_set_domain_attr("i", data_dim=1, type='unstructured' , nvertex=6, i_index=ind_glo) 
     354   CALL xios_set_domain_attr("i",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat) 
     355    
     356   DEALLOCATE(lon, lat, bounds_lon, bounds_lat,ind_glo)  
     357 
     358   dtime%second=1 
     359   CALL xios_set_timestep(dtime) 
     360 
     361   CALL getin('etat0',etat0_type) 
     362   CALL getin('read_metric', read_metric_)  
     363    
     364   CALL xios_set_file_attr('start', enabled=.FALSE.) 
     365   IF (TRIM(etat0_type)=='start_file' .AND. read_metric_) THEN 
     366     CALL xios_set_file_attr('start', enabled=.TRUE.) 
     367   ENDIF 
     368      
     369 
     370   CALL xios_close_context_definition() 
     371!$OMP END MASTER 
     372!$OMP BARRIER 
     373    
     374 END SUBROUTINE xios_init_write_field_input  
    277375   
    278376 SUBROUTINE xios_write_field(name,field) 
    … …  
    9551053  END SUBROUTINE xios_set_context 
    9561054 
     1055 SUBROUTINE xios_set_context_input 
     1056 IMPLICIT NONE    
     1057  TYPE(xios_context) :: ctx_hdl 
     1058 
     1059!$OMP MASTER  
     1060   CALL xios_get_handle("icosagcm_input",ctx_hdl) 
     1061   CALL xios_set_current_context(ctx_hdl) 
     1062!$OMP END MASTER 
     1063 
     1064  END SUBROUTINE xios_set_context_input 
     1065 
    9571066 
    9581067#else 
    … …  
    10201129  SUBROUTINE xios_set_context 
    10211130  END SUBROUTINE xios_set_context 
     1131 
     1132  SUBROUTINE xios_set_context_input 
     1133  END SUBROUTINE xios_set_context_input 
    10221134   
    10231135  SUBROUTINE xios_set_fieldgroup_attr(name,enabled,freq_op) 

  • codes/icosagcm/trunk/src/sphere/geometry.f90

    r814 r882  
    295295    CALL update_circumcenters 
    296296 
    297     DO ind=1,ndomain 
    298       IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE 
    299       d=>domain(ind) 
    300       CALL swap_dimensions(ind) 
    301       CALL swap_geometry(ind) 
    302       DO j=jj_begin,jj_end 
    303         DO i=ii_begin,ii_end 
    304           n=(j-1)*iim+i 
    305           DO k=0,5 
    306             x1(:) = xyz_v(n+z_pos(k+1),:) 
    307             x2(:) = d%vertex(:,k,i,j)  
    308             IF (norm(x1-x2)>1e-10) THEN 
    309               PRINT*,"vertex diff ",ind,i,j,k 
    310               PRINT*,x1 
    311               PRINT*,x2 
    312             ENDIF 
    313           ENDDO 
    314         ENDDO 
    315       ENDDO 
    316     ENDDO 
     297!ym you get discrepency if you try to read metric from file, vertex are not saved in file 
     298!    DO ind=1,ndomain 
     299!      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE 
     300!      d=>domain(ind) 
     301!      CALL swap_dimensions(ind) 
     302!      CALL swap_geometry(ind) 
     303!      DO j=jj_begin,jj_end 
     304!        DO i=ii_begin,ii_end 
     305!          n=(j-1)*iim+i 
     306!          DO k=0,5 
     307!            x1(:) = xyz_v(n+z_pos(k+1),:) 
     308!            x2(:) = d%vertex(:,k,i,j)  
     309!            IF (norm(x1-x2)>1e-10) THEN 
     310!              PRINT*,"vertex diff ",ind,i,j,k 
     311!              PRINT*,x1 
     312!              PRINT*,x2 
     313!            ENDIF 
     314!          ENDDO 
     315!        ENDDO 
     316!      ENDDO 
     317!   ENDDO 
    317318     
    318319     
Note: See TracChangeset for help on using the changeset viewer.