Changeset 856


Ignore:
Timestamp:
05/06/19 01:49:12 (5 years ago)
Author:
dubos
Message:

devel : towards Fortran driver for unstructured/LAM meshes

Location:
codes/icosagcm/devel
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/devel/make_python

    r852 r856  
    2626} 
    2727 
    28 function cmd_extract_dysl() 
    29 { 
    30     XT=$DYNAMICO_ROOT/Python/src/kernels_extracted.jin 
    31     rm -f $XT 
    32     touch $XT 
    33     cd $DYNAMICO_ROOT/src 
    34     for F90 in */*.F90 ; do 
    35         sed -n '/BEGIN_DYSL/,/END_DYSL/{//b;p}' $F90 >> $XT 
    36     done 
    37     cd $DYNAMICO_ROOT 
    38 } 
    39  
    4028function update_kernels() 
    4129{ 
    … …  
    4937{ 
    5038    cmd_clean 
    51     cmd_extract_dysl 
    5239    cd $KERNELS 
    5340    ./codegen hexagonal hex_master unstructured 

  • codes/icosagcm/devel/src/base/dimensions.f90

    r533 r856  
    8080     
    8181    TYPE(t_domain),POINTER :: d 
    82      
     82 
     83    IF(.NOT. swap_needed) RETURN 
     84 
    8385    d=>domain(ind) 
    8486      

  • codes/icosagcm/devel/src/dissip/dissip_gcm.f90

    r845 r856  
    123123    CALL getin("niterdivgrad",niterdivgrad) 
    124124 
    125     CALL dissip_constants 
    126     CALL dissip_profile 
     125    IF(grid_type == grid_ico) THEN 
     126       CALL dissip_constants 
     127       CALL dissip_profile 
     128    END IF 
     129 
    127130  END SUBROUTINE init_dissip 
    128131 

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

    r846 r856  
    4242  !$OMP PARALLEL   
    4343    CALL switch_omp_no_distrib_level 
    44     IF(grid_type == grid_ico) CALL compute_geometry 
     44    IF(grid_type == grid_ico) THEN 
     45       CALL compute_geometry 
     46    ELSE 
     47       swap_needed=.FALSE. ! do this inside an OMP PARALLEL loop to affect all threads 
     48    END IF 
    4549    CALL check_total_area 
    4650   

  • codes/icosagcm/devel/src/initial/etat0_collocated.f90

    r848 r856  
    11MODULE etat0_collocated_mod 
    22  USE icosa 
     3  USE grid_param 
    34  USE omp_para, ONLY : is_omp_level_master 
    45  USE caldyn_vars_mod, ONLY : hydrostatic 
    … …  
    67  PRIVATE 
    78 
    8   LOGICAL :: autoinit_mass, autoinit_NH 
    99  CHARACTER(len=255),SAVE :: etat0_type 
    10 !$OMP THREADPRIVATE(autoinit_mass, autoinit_NH, etat0_type) 
     10!$OMP THREADPRIVATE(etat0_type) 
    1111 
    1212    PUBLIC :: etat0_type, etat0_collocated 
    … …  
    3636    TYPE(t_field),POINTER :: f_q(:) 
    3737   
    38     TYPE(t_field),POINTER,SAVE :: f_temp(:) 
     38    TYPE(t_field),POINTER,SAVE :: f_temp(:) ! SAVE => all threads see the same f_temp 
    3939    REAL(rstd),POINTER :: ps(:) 
    4040    REAL(rstd),POINTER :: mass(:,:) 
    … …  
    6464      q=f_q(ind) 
    6565 
    66       IF( TRIM(etat0_type)=='williamson91.6' ) THEN 
    67          CALL compute_etat0_collocated_hex(ps,mass, phis, theta_rhodz(:,:,1), u, geopot, W, q) 
    68       ELSE 
    69          CALL compute_etat0_collocated_hex(ps,mass, phis, temp, u, geopot, W, q) 
    70       ENDIF 
    71  
    72       IF( TRIM(etat0_type)/='williamson91.6' ) CALL compute_temperature2entropy(ps,temp,q,theta_rhodz, 1) 
    73      
     66      SELECT CASE(grid_type) 
     67      CASE(grid_ico) 
     68         CALL compute_etat0_collocated_hex(phis, ps, mass, theta_rhodz(:,:,1), u, geopot, W, q) 
     69      CASE(grid_unst) 
     70         CALL compute_etat0_collocated_unst(phis, ps, mass, theta_rhodz(:,:,1), u, geopot, W, q) 
     71      CASE DEFAULT 
     72         STOP 'Unexpected value of grid_type encountered in etat0_collocated.' 
     73      END SELECT 
     74 
    7475    ENDDO 
    7576     
    … …  
    7879  END SUBROUTINE etat0_collocated 
    7980 
    80   SUBROUTINE compute_temperature2entropy(ps,temp,q,theta_rhodz,offset) 
    81     USE icosa 
    82     USE pression_mod 
    83     USE exner_mod 
    84     USE omp_para 
    85     REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    86     REAL(rstd),INTENT(IN)  :: temp(iim*jjm,llm) 
    87     REAL(rstd),INTENT(IN)  :: q(iim*jjm,llm,nqtot) 
    88     REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 
    89     INTEGER,INTENT(IN) :: offset 
    90  
    91     REAL(rstd) :: p(iim*jjm,llm+1) 
    92     REAL(rstd) :: cppd,Rd, mass, p_ij, q_ij,r_ij, chi,nu, entropy, theta 
    93     INTEGER :: i,j,ij,l 
    94  
    95     cppd=cpp 
    96     Rd=kappa*cppd 
    97  
    98     CALL compute_pression(ps,p,offset) 
    99     ! flush p 
    100     DO    l    = ll_begin, ll_end 
    101        DO j=jj_begin-offset,jj_end+offset 
    102           DO i=ii_begin-offset,ii_end+offset 
    103              ij=(j-1)*iim+i 
    104              mass = (p(ij,l)-p(ij,l+1))/g ! dry+moist mass 
    105              p_ij = .5*(p(ij,l)+p(ij,l+1))  ! pressure at full level 
    106              SELECT CASE(caldyn_thermo) 
    107              CASE(thermo_theta) 
    108                 theta = temp(ij,l)*(p_ij/preff)**(-kappa)  
    109                 theta_rhodz(ij,l) = mass * theta 
    110              CASE(thermo_entropy) 
    111                 nu = log(p_ij/preff) 
    112                 chi = log(temp(ij,l)/Treff) 
    113                 entropy = cppd*chi-Rd*nu 
    114                 theta_rhodz(ij,l) = mass * entropy 
    115 !             CASE(thermo_moist) 
    116 !                q_ij=q(ij,l,1) 
    117 !                r_ij=1.-q_ij 
    118 !                mass=mass*(1-q_ij) ! dry mass 
    119 !                nu = log(p_ij/preff) 
    120 !                chi = log(temp(ij,l)/Treff) 
    121 !                entropy = r_ij*(cppd*chi-Rd*nu) + q_ij*(cppv*chi-Rv*nu) 
    122 !                theta_rhodz(ij,l) = mass * entropy                 
    123                 CASE DEFAULT 
    124                    STOP 
    125              END SELECT 
    126           ENDDO 
    127        ENDDO 
    128     ENDDO 
    129   END SUBROUTINE compute_temperature2entropy 
    130  
    131   SUBROUTINE compute_etat0_collocated_hex(ps,mass,phis,temp_i,u, geopot,W, q) 
    132     USE wind_mod 
    133     USE disvert_mod 
     81  SUBROUTINE compute_etat0_collocated_hex(phis,ps,mass,theta_rhodz,u, geopot,W, q) 
    13482    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) 
    13583    REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm) 
     84    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 
    13685    REAL(rstd),INTENT(OUT) :: phis(iim*jjm) 
    137     REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm) 
    13886    REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) 
    13987    REAL(rstd),INTENT(OUT) :: W(iim*jjm,llm+1) 
    … …  
    14189    REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) 
    14290 
    143     REAL(rstd) :: ulon_i(iim*jjm,llm) 
    144     REAL(rstd) :: ulat_i(iim*jjm,llm) 
    145  
    14691    REAL(rstd) :: ps_e(3*iim*jjm) 
     92    REAL(rstd) :: phis_e(3*iim*jjm) 
    14793    REAL(rstd) :: mass_e(3*iim*jjm,llm) 
    148     REAL(rstd) :: phis_e(3*iim*jjm) 
    149     REAL(rstd) :: temp_e(3*iim*jjm,llm) 
    15094    REAL(rstd) :: geopot_e(3*iim*jjm,llm+1) 
    151     REAL(rstd) :: ulon_e(3*iim*jjm,llm) 
    152     REAL(rstd) :: ulat_e(3*iim*jjm,llm) 
    15395    REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot) 
    15496 
    155     INTEGER :: l,i,j,ij 
    156     REAL :: p_ik, v_ik, mass_ik 
    157  
    158     ! For NH geopotential and vertical momentum must be initialized. 
    159     ! Unless autoinit_NH is set to .FALSE. , they will be initialized 
    160     ! to hydrostatic geopotential and zero 
    161     autoinit_mass = .TRUE. 
    162     autoinit_NH = .NOT. hydrostatic 
    16397    w(:,:) = 0 
    164  
    165     CALL compute_etat0_collocated(iim*jjm  , lon_i, lat_i, phis,   ps,   mass,   temp_i, ulon_i, ulat_i, geopot,   q) 
    166     CALL compute_etat0_collocated(3*iim*jjm, lon_e, lat_e, phis_e, ps_e, mass_e, temp_e, ulon_e, ulat_e, geopot_e, q_e) 
    167  
    168     IF(autoinit_mass) CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps 
    169     IF(autoinit_NH) THEN 
    170        geopot(:,1) = phis(:) ! surface geopotential 
    171        DO l = 1, llm 
    172           DO ij=1,iim*jjm 
    173              ! hybrid pressure coordinate 
    174              p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) 
    175              mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g 
    176              ! v=R.T/p, R=kappa*cpp 
    177              v_ik = kappa*cpp*temp_i(ij,l)/p_ik 
    178              geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g 
    179           END DO 
    180        END DO 
    181     END IF 
    182  
    183     CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u) 
     98    CALL compute_etat0_collocated(iim*jjm  , lon_i, lat_i, phis,   ps,   mass,   theta_rhodz, geopot,   q) 
     99    CALL compute_etat0_collocated(3*iim*jjm, lon_e, lat_e, phis_e, ps_e, mass_e, mass_e, geopot_e, q_e, ep_e, u) 
    184100 
    185101  END SUBROUTINE compute_etat0_collocated_hex 
    186102 
    187   SUBROUTINE compute_etat0_collocated(ngrid, lon, lat, phis, ps, mass, temp, ulon, ulat, geopot, q) 
     103  SUBROUTINE compute_etat0_collocated_unst(phis,ps,mass,theta_rhodz,u, geopot,W, q) 
     104    REAL(rstd),INTENT(INOUT) :: ps(primal_num) 
     105    REAL(rstd),INTENT(INOUT) :: mass(llm,primal_num) 
     106    REAL(rstd),INTENT(OUT) :: theta_rhodz(llm,primal_num) 
     107    REAL(rstd),INTENT(OUT) :: phis(primal_num) 
     108    REAL(rstd),INTENT(OUT) :: u(llm, edge_num) 
     109    REAL(rstd),INTENT(OUT) :: W(llm+1, primal_num) 
     110    REAL(rstd),INTENT(OUT) :: geopot(llm+1, primal_num) 
     111    REAL(rstd),INTENT(OUT) :: q(llm, primal_num, nqtot) 
     112    ! The 2D/3D arrays below have the shape expected by compute_etat0_collocated 
     113    REAL(rstd) :: ps_e(edge_num) 
     114    REAL(rstd) :: phis_e(edge_num) 
     115    REAL(rstd) :: u_e(edge_num, llm) 
     116    REAL(rstd) :: mass_i(primal_num, llm), theta_rhodz_i(primal_num, llm), mass_e(edge_num, llm) 
     117    REAL(rstd) :: geopot_i(edge_num, llm+1), geopot_e(edge_num, llm+1) 
     118    REAL(rstd) :: q_i(primal_num, llm, nqtot), q_e(edge_num, llm, nqtot) 
     119    INTEGER :: iq 
     120 
     121    w(:,:) = 0 
     122    CALL compute_etat0_collocated(primal_num  , lon_i, lat_i, phis,   ps,   mass_i,   theta_rhodz_i, geopot_i,   q_i) 
     123    CALL compute_etat0_collocated(edge_num, lon_e, lat_e, phis_e, ps_e, mass_e, mass_e, geopot_e, q_e, ep_e, u_e) 
     124    mass = TRANSPOSE(mass_i) 
     125    theta_rhodz = TRANSPOSE(theta_rhodz_i) 
     126    geopot = TRANSPOSE(geopot_i) 
     127    u = TRANSPOSE(u_e) 
     128    DO iq=1,nqtot 
     129       q(:,:,iq) = TRANSPOSE(q_i(:,:,iq)) 
     130    END DO 
     131  END SUBROUTINE compute_etat0_collocated_unst 
     132   
     133  SUBROUTINE compute_etat0_collocated(ngrid, lon, lat, phis, ps, mass, theta_rhodz, geopot, q, ep, uperp) 
    188134    USE etat0_isothermal_mod, ONLY : compute_isothermal => compute_etat0 
    189135    USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0 
    … …  
    199145    USE etat0_dcmip2016_cyclone_mod, ONLY : compute_dcmip2016_cyclone => compute_etat0 
    200146    USE etat0_dcmip2016_supercell_mod, ONLY : compute_dcmip2016_supercell => compute_etat0 
     147    USE disvert_mod, ONLY : ptop, ap, bp, mass_ak, mass_bk, mass_dak, mass_dbk 
    201148    INTEGER :: ngrid 
    202149    REAL(rstd),INTENT(IN)    :: lon(ngrid), lat(ngrid) 
    203150    REAL(rstd),INTENT(INOUT) :: ps(ngrid) 
    204151    REAL(rstd),INTENT(INOUT) :: mass(ngrid,llm) 
     152    REAL(rstd),INTENT(OUT)   :: theta_rhodz(ngrid,llm) 
    205153    REAL(rstd),INTENT(OUT)   :: phis(ngrid) 
    206     REAL(rstd),INTENT(OUT)   :: temp(ngrid,llm) 
    207     REAL(rstd),INTENT(OUT)   :: ulon(ngrid,llm) 
    208     REAL(rstd),INTENT(OUT)   :: ulat(ngrid,llm) 
    209154    REAL(rstd),INTENT(OUT)   :: geopot(ngrid,llm+1) 
    210155    REAL(rstd),INTENT(OUT)   :: q(ngrid,llm,nqtot) 
     156    REAL(rstd),INTENT(OUT), OPTIONAL   :: ep(ngrid,3), uperp(ngrid,llm) 
     157 
     158    REAL(rstd)   :: ulon(ngrid,llm) 
     159    REAL(rstd)   :: ulat(ngrid,llm) 
     160    REAL(rstd)   :: temp(ngrid,llm) 
     161 
     162    LOGICAL :: autoinit_mass, autoinit_NH 
     163    INTEGER :: ij,l 
     164    REAL(rstd) :: clat, slat, clon, slon, u3d(3), p_ik, mass_ik, v_ik 
     165    REAL(rstd) :: q_ik, r_ik, chi, entropy, theta, log_p_preff 
     166 
     167    ! For NH, geopotential and vertical momentum must be initialized. 
     168    ! Unless autoinit_NH is set to .FALSE. , they will be initialized 
     169    ! to hydrostatic geopotential and zero 
     170    autoinit_mass = .TRUE. 
     171    autoinit_NH = .NOT. hydrostatic 
    211172 
    212173    SELECT CASE (TRIM(etat0_type)) 
    … …  
    232193!       autoinit_NH = .FALSE. ! compute_bubble initializes geopot 
    233194    CASE('williamson91.6') 
    234        CALL compute_w91_6(ngrid, lon, lat, phis, mass(:,1), temp(:,1), ulon(:,1), ulat(:,1)) 
     195       CALL compute_w91_6(ngrid, lon, lat, phis, mass(:,1), theta_rhodz(:,1), ulon(:,1), ulat(:,1)) 
    235196       autoinit_mass = .FALSE. ! do not overwrite mass 
    236197    CASE('dcmip2016_baroclinic_wave') 
    … …  
    242203    END SELECT 
    243204 
     205    IF(PRESENT(uperp)) THEN 
     206       DO l = 1, llm 
     207          DO ij=1,ngrid 
     208             clat = COS(lat(ij)) 
     209             slat = SIN(lat(ij)) 
     210             clon = COS(lon(ij)) 
     211             slon = SIN(lon(ij)) 
     212             u3d(1) = -clon*slat*ulat(ij,l) - slon*ulon(ij,l) 
     213             u3d(2) = -slon*slat*ulat(ij,l) + clon*ulon(ij,l) 
     214             u3d(3) =       clat*ulat(ij,l)     
     215             uperp(ij,l) = SUM(u3d*ep(ij,:)) 
     216          END DO 
     217       END DO 
     218    END IF 
     219 
     220    IF(autoinit_mass) THEN 
     221       DO l = 1, llm 
     222          DO ij=1,ngrid 
     223             mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g 
     224             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) 
     225             mass(ij,l) = mass_ik 
     226             SELECT CASE(caldyn_thermo) 
     227             CASE(thermo_theta) 
     228                theta = temp(ij,l)*(p_ik/preff)**(-kappa) 
     229                theta_rhodz(ij,l) = mass_ik * theta 
     230             CASE(thermo_entropy) 
     231                log_p_preff = log(p_ik/preff) 
     232                chi = log(temp(ij,l)/Treff) 
     233                entropy = cpp*chi-Rd*log_p_preff 
     234                theta_rhodz(ij,l) = mass_ik * entropy 
     235             CASE(thermo_moist) 
     236                q_ik=q(ij,l,1) 
     237                r_ik=1.-q_ik 
     238                mass_ik = mass_ik*(1.-q_ik) ! dry mass 
     239                log_p_preff = log(p_ik/preff) 
     240                chi = log(temp(ij,l)/Treff) 
     241                entropy = r_ik*(cpp*chi-Rd*nu) + q_ik*(cppv*chi-Rv*log_p_preff) 
     242                theta_rhodz(ij,l) = mass_ik * entropy 
     243             CASE DEFAULT 
     244                STOP 
     245             END SELECT 
     246          END DO 
     247       END DO 
     248    END IF 
     249 
     250    IF(autoinit_NH) THEN 
     251       geopot(:,1) = phis(:) ! surface geopotential 
     252       DO l = 1, llm 
     253          DO ij=1,ngrid 
     254             ! hybrid pressure coordinate 
     255             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) 
     256             mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g 
     257             ! v=R.T/p, R=kappa*cpp 
     258             v_ik = kappa*cpp*temp(ij,l)/p_ik 
     259             geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g 
     260          END DO 
     261       END DO 
     262    END IF 
     263 
    244264  END SUBROUTINE compute_etat0_collocated 
    245265 

  • codes/icosagcm/devel/src/output/output_field.f90

    r714 r856  
    2020  SUBROUTINE output_field_init 
    2121  USE getin_mod 
     22  USE grid_param 
    2223  IMPLICIT NONE 
    2324 
    … …  
    3536     
    3637    IF (xios_output) THEN 
    37       CALL xios_init_write_field 
     38       IF(grid_type == grid_ico) CALL xios_init_write_field 
    3839    ENDIF 
    3940  END SUBROUTINE output_field_init 

  • codes/icosagcm/devel/src/parallel/domain.f90

    r821 r856  
    6565  INTEGER,SAVE :: ndomain 
    6666  INTEGER,SAVE :: ndomain_glo 
     67 
     68  LOGICAL :: swap_needed=.TRUE. ! .FALSE. if a thread always computes on the same domain 
     69  !$OMP THREADPRIVATE(swap_needed) 
    6770 
    6871  TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain(:) 

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

    r602 r856  
    142142  SUBROUTINE swap_geometry(ind) 
    143143  USE field_mod 
     144  USE domain_mod, ONLY : swap_needed 
    144145  IMPLICIT NONE 
    145146    INTEGER,INTENT(IN) :: ind 
     147 
     148    IF(.NOT. swap_needed) RETURN 
     149 
    146150!!$OMP MASTER 
    147151    Ai=geom%Ai(ind) 

  • codes/icosagcm/devel/src/transport/advect_tracer.f90

    r592 r856  
    2727  SUBROUTINE init_advect_tracer 
    2828    USE omp_para 
     29    USE grid_param 
    2930    REAL(rstd),POINTER :: tangent(:,:) 
    3031    REAL(rstd),POINTER :: normal(:,:) 
    … …  
    4243    CALL allocate_field(f_wq, field_t, type_real, llm+1, name='wq') 
    4344     
    44     DO ind=1,ndomain 
    45        IF (.NOT. assigned_domain(ind)) CYCLE 
    46        CALL swap_dimensions(ind) 
    47        CALL swap_geometry(ind) 
    48        normal=f_normal(ind) 
    49        tangent=f_tangent(ind) 
    50        sqrt_leng=f_sqrt_leng(ind) 
    51        IF (is_omp_level_master) CALL init_advect(normal,tangent,sqrt_leng) 
    52     END DO 
     45    IF(grid_type == grid_ico) THEN 
     46       DO ind=1,ndomain 
     47          IF (.NOT. assigned_domain(ind)) CYCLE 
     48          CALL swap_dimensions(ind) 
     49          CALL swap_geometry(ind) 
     50          normal=f_normal(ind) 
     51          tangent=f_tangent(ind) 
     52          sqrt_leng=f_sqrt_leng(ind) 
     53          IF (is_omp_level_master) CALL init_advect(normal,tangent,sqrt_leng) 
     54       END DO 
     55    END IF 
    5356 
    5457  END SUBROUTINE init_advect_tracer 

  • codes/icosagcm/devel/src/unstructured/init_unstructured.f90

    r830 r856  
    9696    USE ioipsl 
    9797    USE field_mod 
     98    USE geometry, ONLY : lon_i, lat_i, lon_e, lat_e, ep_e 
    9899    IMPLICIT NONE 
    99100 
    100101    !!-------------Declare local variables------------------- 
    101     CHARACTER(LEN=*),PARAMETER :: f="mesh_information.nc" 
    102     INTEGER :: id_nc, ierr, status, descriptionLength 
     102    CHARACTER(LEN=*),PARAMETER :: f="input/mesh_information.nc" 
     103    INTEGER :: id_nc, ierr, status, descriptionLength, ij 
    103104    CHARACTER(LEN= 80) :: description 
    104  
     105    REAL(rstd), ALLOCATABLE :: angle_e(:) 
     106    REAL(rstd) :: clon, slon, clat, slat, & ! COS/SIN of lon/lat 
     107         elon(3), elat(3) ! lon/lat unit vectors 
    105108    print *,"------------------ READING FILE " , f, "----------------------- " 
    106109    !open and read the input file 
    … …  
    167170    ALLOCATE(dual_vertex, source = Idata_read2) 
    168171 
     172    CALL read_from_gridfile(id_nc, 'float', 'lon_i') 
     173    ALLOCATE(lon_i, source = Ddata_read1) 
     174    CALL read_from_gridfile(id_nc, 'float', 'lat_i') 
     175    ALLOCATE(lat_i, source = Ddata_read1) 
     176    CALL read_from_gridfile(id_nc, 'float', 'lon_e') 
     177    ALLOCATE(lon_e, source = Ddata_read1) 
     178    CALL read_from_gridfile(id_nc, 'float', 'lat_e') 
     179    ALLOCATE(lat_e, source = Ddata_read1) 
     180    CALL read_from_gridfile(id_nc, 'float', 'angle_e') 
     181    ALLOCATE(angle_e, source = Ddata_read1) 
     182 
    169183    CALL free_data_read ! free buffers after reading all data from grid file 
    170184 
    … …  
    176190    max_trisk_deg = SIZE(trisk,1) 
    177191 
     192    ! now post-process some of the data we just read 
     193    ALLOCATE(ep_e(edge_num,3)) 
     194    DO ij=1,edge_num 
     195       clon = COS(lon_e(ij)) 
     196       slon = SIN(lon_e(ij)) 
     197       clat = COS(lat_e(ij)) 
     198       slat = SIN(lat_e(ij)) 
     199       ! x,y,z = clat*clon, clat*slon, slat 
     200       elon = (/ -slon, clon, 0. /) 
     201       elat = (/ -slat*clon, -slat*slon, clat /) 
     202       ep_e(ij,:) = COS(angle_e(ij))*elon + SIN(angle_e(ij))*elat 
     203    END DO 
     204     
     205    DEALLOCATE(angle_e) 
    178206  END SUBROUTINE read_dump_partition 
    179207 
Note: See TracChangeset for help on using the changeset viewer.