Changes in / [20:30]


Ignore:
Location:
/codes/icosagcm/trunk
Files:
12 added
1 deleted
25 edited

Legend:

Unmodified
Added
Removed

  • /codes/icosagcm/trunk/arch/arch-X64_TITANE.fcm

    r20 r30  
    44%MAKE                gmake 
    55%FPP_FLAGS           -P -traditional 
    6 %FPP_DEF             NC_DOUBLE BLAS SGEMV=DGEMV SGEMM=DGEMM FFT_MKL 
    7 %BASE_FFLAGS        -openmp  -i4 -r8 -automatic -align all -I${MKLROOT}/include -I$NETCDF_INC_DIR 
    8 %PROD_FFLAGS         -O3 
     6%FPP_DEF             KEY_NONE  
     7%BASE_FFLAGS         -i4 -r8 -automatic -align all -I${MKLROOT}/include 
     8%PROD_FFLAGS         -O3 -openmp 
    99%DEV_FFLAGS          -p -g -O3 -traceback -fp-stack-check -ftrapuv 
    10 %DEBUG_FFLAGS        -p -g -traceback 
     10%DEBUG_FFLAGS        -p -g -traceback -check bounds 
    1111%MPI_FFLAGS 
    1212%OMP_FFLAGS          -openmp 
    13 %BASE_LD             -openmp -i4 -r8 -automatic $MKL_LIBS -L$NETCDF_LIB_DIR -lnetcdf -lnetcdff 
     13%BASE_LD             -openmp -i4 -r8 -automatic $MKL_LIBS 
    1414%MPI_LD 
    1515%OMP_LD              -openmp 

  • /codes/icosagcm/trunk/arch/arch-X64_TITANE.path

    r20 r30  
    1 NETCDF_LIBDIR="$NETCDF_LIB_DIR -lnetcdff" 
    2 NETCDF_INCDIR=$NETCDF_INC_DIR 
    3 IOIPSL_INCDIR=$LMDGCM/../../lib 
    4 IOIPSL_LIBDIR=$LMDGCM/../../lib 
    5 ORCH_INCDIR=$LMDGCM/../../lib 
    6 ORCH_LIBDIR=$LMDGCM/../../lib 
    7 OASIS_INCDIR=$LMDGCM/../../prism/X64/build/lib/psmile.$couple 
    8 OASIS_LIBDIR=$LMDGCM/../../prism/X64/lib 
    9 INCA_LIBDIR=$LMDGCM/../INCA3/config/lib 
    10 INCA_INCDIR=$LMDGCM/../INCA3/config/lib 
     1NETCDF_INCDIR="-I $NETCDF_INC_DIR" 
     2NETCDF_LIBDIR="-L $NETCDF_LIB_DIR" 
     3NETCDF_LIB="-lnetcdff -lnetcdf" 
    114 
     5MPI_INCDIR="" 
     6MPI_LIBDIR="" 
     7MPI_LIB="" 
     8 
     9HDF5_INCDIR="-I $HDF5_INC_DIR" 
     10HDF5_LIBDIR="-L $HDF5_LIB_DIR" 
     11HDF5_LIB="-lhdf5_hl -lhdf5 -lhdf5 -lz -lcurl" 

  • /codes/icosagcm/trunk/bld.cfg

    r20 r30  
    2626bld::tool::fc        %COMPILER 
    2727bld::tool::ld        %LINK 
    28 bld::tool::ldflags   %LD_FLAGS  
     28bld::tool::ldflags   %LD_FLAGS %LIB  
    2929bld::tool::fflags    %FFLAGS  
     30bld::tool::fppkeys   %CPP_KEY %FPP_DEF 
     31 
    3032# Pre-process code before analysing dependencies 
    31 bld::pp              false 
     33bld::pp              true 
    3234 
    3335bld::excl_dep        use::netcdf 
    3436bld::excl_dep        use::omp_lib 
     37bld::excl_dep        inc::mpif.h 
     38 
    3539bld::tool::SHELL   /bin/bash 
    3640 

  • /codes/icosagcm/trunk/make_icosa

    r20 r30  
    77 
    88arch_defined="FALSE" 
     9parallel_defined="FALSE" 
    910arch="" 
     11parallel="none" 
     12CPP_KEY="CPP_NONE"  
    1013 
    1114while (($# > 0)) 
    … …  
    3336          arch=$2 ; arch_defined="TRUE"; shift ; shift ;; 
    3437 
     38      "-parallel") 
     39          parallel=$2 ; parallel_defined="TRUE"; shift ; shift ;; 
     40 
    3541      *) 
    3642          code="$1" ; shift ;; 
    … …  
    3844done 
    3945 
     46rm -f .void_file 
     47echo > .void_file 
     48rm -rf .void_dir 
     49mkdir .void_dir 
     50 
    4051if [[ "$arch_defined" == "TRUE" ]] 
    4152then 
    4253  rm -f arch.path 
    4354  rm -f arch.fcm 
     55  rm -f arch.env 
    4456  ln -s arch/arch-${arch}.path ./arch.path 
    4557  ln -s arch/arch-${arch}.fcm  ./arch.fcm 
     58  if test -f arch/arch-${arch}.env 
     59    then 
     60      ln -s arch/arch-${arch}.env arch.env 
     61    else 
     62      ln -s .void_file arch.env 
     63    fi 
     64  source arch.env 
    4665  source arch.path 
    4766else 
    … …  
    6180fi 
    6281 
     82if [[ "$parallel" == "mpi" ]] 
     83then 
     84  CPP_KEY="$CPP_KEY CPP_USING_MPI" 
     85elif [[ "$parallel" == "none" ]] 
     86then 
     87  parallel="none" 
     88else  
     89  echo "-parallel value $parallel is invalid, only permited <none> or <mpi>" 
     90  exit 1 
     91fi 
     92 
     93ICOSA_LIB="$NETCDF_LIBDIR $HDF5_LIBDIR $NETCDF_LIB $HDF5_LIB" 
     94 
    6395rm -f config.fcm 
    6496 
    65 echo "%COMPIL_FFLAGS $COMPIL_FFLAGS" >> config.fcm 
     97echo "%COMPIL_FFLAGS $COMPIL_FFLAGS $NETCDF_INCDIR" >> config.fcm 
     98echo "%CPP_KEY $CPP_KEY" >> config.fcm 
     99echo "%LIB $ICOSA_LIB">> config.fcm 
    66100 
    67101./build 

  • /codes/icosagcm/trunk/run_adv.def

    r20 r30  
     1#------------------------------- Mesh --------------------------------- 
     2 
    13# Number of subdivision on a main triangle (nbp) : integer (default=40) 
    24nbp=20 
    35 
    4 # Number of vertical layer (llm) : integer (default=19) 
     6# optim_it : mesh optimisation : number of iteration : integer (default=0) 
     7optim_it=0 
     8 
     9# Number of vertical layers (llm) : integer (default=19) 
    510llm=19 
     11 
     12# disvert : vertical discretisation : string (default='std') : std, ncar 
     13disvert=ncar 
     14 
     15#---------------------------------- Misc -------------------------------- 
    616 
    717# number of tracer (nqtot) : integer (default 1) 
    … …  
    1121dt = 180. 
    1222 
    13 # scheme type : string ( default='adam_bashforth') euler, leapfrog, leapfrog_matsuno, adam_bashforth) 
     23# number of timestep (default 100) 
     24itaumax = 38400 
     25 
     26# output field period : integer (default none) 
     27write_period=3600 
     28 
     29# etat0 : initial state : string (default=jablonowsky06) :  
     30# jablonowsky06, academic, ncar 
     31etat0=ncar 
     32 
     33# caldyn : computation type for gcm equation : string (default=gcm) : gcm, adv 
     34caldyn=adv 
     35 
     36#------------------------------ Dynamics -------------------------------- 
     37 
     38# scheme type : string ( default='adam_bashforth') 
     39# euler, leapfrog, leapfrog_matsuno, adam_bashforth) 
    1440scheme = euler 
    1541 
    1642# matsuno period : integer ( default=5) 
    1743matsuno_period = 5 
    18  
    19 # number of timestep (default 100) 
    20 itaumax = 38400 
    2144 
    2245# dissipation time graddiv : real (default=5000) 
    … …  
    2548 
    2649# dissipation time nxgradrot (default=5000) 
    27 tetagrot = 1800 
    28 tau_gradrot = 1800 
    29  
     50tetagrot=1800 
     51tau_gradrot=1800 
    3052tau_divgrad=1800 
    31  
    32 # output field period : integer (default none) 
    33 write_period=3600 
    34  
    35  
    36 #NCAR testcase : integer ( default=5) 
    37 ncar_test_case=1 
    38  
    39 # disvert : vertical discretisation : string (default='std') : std, ncar 
    40 disvert=ncar 
    41  
    42 # optim_it : mesh optimisation : number of iteration : integer (default=0) 
    43 optim_it=0 
    44  
    45 # initial_state 
    46  
    4753 
    4854# guided_type : string (default=none) : none, ncar 
    4955guided_type=ncar 
    5056 
    51 # etat0 : initial state : string (default=jablonowsky06) : jablonowsky06, academic, ncar 
    52 etat0=ncar 
     57#-------------------- parameters for NCAR test cases ------------------------ 
    5358 
    54 # caldyn : computation type for gcm equation : string (default=gcm) : gcm, adv 
    55 caldyn=adv 
     59# NCAR advection test, initial tracer : string ( default='cos_bell') 
     60# const, slotted_cyl, cos_bell, dbl_cos_bell_q1, dbl_cos_bell_q2, complement, hadley 
     61ncar_adv_shape=complement 
     62 
     63# NCAR advection test, wind field : string (default='deform') : solid, deform, hadley 
     64ncar_adv_wind=deform 
     65 
     66# ncar_dz : model layer thickness in meters: real (default=400)  
     67# used if disvert=ncar 
     68ncar_dz=400 
     69 
     70# ncar_T0 : reference temperature for NCAR test cases : real (default=300)  
     71# also used by disvert if disvert=ncar 
     72ncar_T0=300 
     73 
     74# ncar_p0 : reference pressure for NCAR test cases : real (default=1e5)  
     75# also used by disvert if disvert=ncar 
     76ncar_p0=1e5 
     77 
     78# ncar_disvert_c : exponent for B(eta) : integer (default=1) 
     79# used by disvert if disvert=ncar 
     80ncar_disvert_c=1 

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

    r20 r30  
    1       MODULE advect_mod 
    2    USE icosa 
    3         IMPLICIT NONE  
    4  
    5  
    6   TYPE(t_field),POINTER :: f_normal(:) 
    7   TYPE(t_field),POINTER :: f_tangent(:) 
    8   TYPE(t_field),POINTER ::  f_gradq3d(:) 
    9   REAL(rstd),POINTER :: gradq3d(:,:,:)  
    10   REAL(rstd),POINTER :: normal(:,:) 
    11   REAL(rstd),POINTER :: tangent(:,:) 
    12                          
     1MODULE advect_mod 
     2  USE icosa 
     3  IMPLICIT NONE  
     4 
    135CONTAINS 
    14           
    15          
    16   SUBROUTINE allocate_advect 
    17   USE field_mod 
    18   USE domain_mod 
    19   USE dimensions 
    20   USE geometry 
    21   USE metric 
    22   IMPLICIT NONE 
    23  
    24     CALL allocate_field(f_normal,field_u,type_real,3) 
    25     CALL allocate_field(f_tangent,field_u,type_real,3) 
    26     CALL allocate_field(f_gradq3d,field_t,type_real,llm,3) 
    27  
    28   END SUBROUTINE allocate_advect 
    29 !========================================================================== 
    30   SUBROUTINE swap_advect(ind) 
    31   USE domain_mod 
    32   USE dimensions 
    33   USE geometry 
    34   USE metric 
    35   IMPLICIT NONE 
    36   INTEGER,INTENT(IN) :: ind 
    37    
    38       normal=f_normal(ind) 
    39       tangent=f_tangent(ind) 
    40       gradq3d = f_gradq3d(ind)  
    41  
    42   END SUBROUTINE swap_advect 
    43 !========================================================================== 
    44      
    45   SUBROUTINE init_advect 
    46   USE domain_mod 
    47   USE dimensions 
    48   USE geometry 
    49   USE metric 
    50   USE vector 
    51   IMPLICIT NONE 
    52    INTEGER :: ind,i,j,n 
    53     
    54    CALL allocate_advect 
    55     
     6 
     7  !========================================================================== 
     8 
     9  SUBROUTINE init_advect(normal,tangent) 
     10    USE domain_mod 
     11    USE dimensions 
     12    USE geometry 
     13    USE metric 
     14    USE vector 
     15    IMPLICIT NONE 
     16    REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3) 
     17    REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3) 
     18 
     19    INTEGER :: i,j,n 
     20 
    5621    DO j=jj_begin-1,jj_end+1 
    57       DO i=ii_begin-1,ii_end+1 
    58         n=(j-1)*iim+i 
    59              
    60         CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),normal(n+u_right,:)) 
    61         normal(n+u_right,:)=normal(n+u_right,:)/sqrt(sum(normal(n+u_right,:)**2)+1e-50)*ne(n,right) 
    62           
    63         CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),normal(n+u_lup,:)) 
    64          normal(n+u_lup,:)=normal(n+u_lup,:)/sqrt(sum(normal(n+u_lup,:)**2)+1e-50)*ne(n,lup) 
    65          
    66         CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),normal(n+u_ldown,:))     
    67         normal(n+u_ldown,:)=normal(n+u_ldown,:)/sqrt(sum(normal(n+u_ldown,:)**2)+1e-50)*ne(n,ldown) 
    68                        
    69         tangent(n+u_right,:)=xyz_v(n+z_rup,:)-xyz_v(n+z_rdown,:)  
    70         tangent(n+u_right,:)=tangent(n+u_right,:)/sqrt(sum(tangent(n+u_right,:)**2)+1e-50) 
    71               
    72         tangent(n+u_lup,:)=xyz_v(n+z_lup,:)-xyz_v(n+z_up,:)  
    73         tangent(n+u_lup,:)=tangent(n+u_lup,:)/sqrt(sum(tangent(n+u_lup,:)**2)+1e-50) 
    74       
    75         tangent(n+u_ldown,:)=xyz_v(n+z_down,:)-xyz_v(n+z_ldown,:) 
    76         tangent(n+u_ldown,:)=tangent(n+u_ldown,:)/sqrt(sum(tangent(n+u_ldown,:)**2)+1e-50) 
    77       END DO 
    78      ENDDO       
     22       DO i=ii_begin-1,ii_end+1 
     23          n=(j-1)*iim+i 
     24 
     25          CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),normal(n+u_right,:)) 
     26          normal(n+u_right,:)=normal(n+u_right,:)/sqrt(sum(normal(n+u_right,:)**2)+1e-50)*ne(n,right) 
     27 
     28          CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),normal(n+u_lup,:)) 
     29          normal(n+u_lup,:)=normal(n+u_lup,:)/sqrt(sum(normal(n+u_lup,:)**2)+1e-50)*ne(n,lup) 
     30 
     31          CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),normal(n+u_ldown,:))   
     32          normal(n+u_ldown,:)=normal(n+u_ldown,:)/sqrt(sum(normal(n+u_ldown,:)**2)+1e-50)*ne(n,ldown) 
     33 
     34          tangent(n+u_right,:)=xyz_v(n+z_rup,:)-xyz_v(n+z_rdown,:)  
     35          tangent(n+u_right,:)=tangent(n+u_right,:)/sqrt(sum(tangent(n+u_right,:)**2)+1e-50) 
     36 
     37          tangent(n+u_lup,:)=xyz_v(n+z_lup,:)-xyz_v(n+z_up,:)  
     38          tangent(n+u_lup,:)=tangent(n+u_lup,:)/sqrt(sum(tangent(n+u_lup,:)**2)+1e-50) 
     39 
     40          tangent(n+u_ldown,:)=xyz_v(n+z_down,:)-xyz_v(n+z_ldown,:) 
     41          tangent(n+u_ldown,:)=tangent(n+u_ldown,:)/sqrt(sum(tangent(n+u_ldown,:)**2)+1e-50) 
     42       END DO 
     43    ENDDO 
     44 
    7945  END SUBROUTINE init_advect 
    80 !======================================================================================= 
    81  
    82  
    83  
    84 !======================================================================================= 
    85   SUBROUTINE advect1(qi) 
    86   USE domain_mod 
    87   USE dimensions 
    88   USE geometry 
    89   USE metric 
    90   IMPLICIT NONE 
    91     REAL(rstd),INTENT(IN) :: qi(iim*jjm,llm) 
     46 
     47  !======================================================================================= 
     48 
     49  SUBROUTINE compute_gradq3d(qi,gradq3d) 
     50    USE domain_mod 
     51    USE dimensions 
     52    USE geometry 
     53    USE metric 
     54    IMPLICIT NONE 
     55    REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm) 
     56    REAL(rstd),INTENT(OUT) :: gradq3d(iim*jjm,llm,3)  
    9257    REAL(rstd) :: maxq,minq,minq_c,maxq_c  
    9358    REAL(rstd) :: alphamx,alphami,alpha ,maggrd,leng 
    … …  
    9762    REAL(rstd) :: ar 
    9863    INTEGER :: i,j,n,k,ind,l 
    99 !========================================================================================== GRADIENT  
    100        Do l = 1,llm  
    101         DO j=jj_begin-1,jj_end+1 
    102          DO i=ii_begin-1,ii_end+1 
    103           n=(j-1)*iim+i    
    104           CALL gradq(n,n+t_rup,n+t_lup,n+z_up,qi(:,l),gradtri(n+z_up,l,:),arr(n+z_up)) 
    105           CALL gradq(n,n+t_ldown,n+t_rdown,n+z_down,qi(:,l),gradtri(n+z_down,l,:),arr(n+z_down)) 
    106          END DO 
    107         END DO 
    108        END DO 
    109  
    110 !      Do l =1,llm 
    111         DO j=jj_begin,jj_end 
     64    !========================================================================================== GRADIENT  
     65    Do l = 1,llm  
     66       DO j=jj_begin-1,jj_end+1 
     67          DO i=ii_begin-1,ii_end+1 
     68             n=(j-1)*iim+i    
     69             CALL gradq(n,n+t_rup,n+t_lup,n+z_up,qi(:,l),gradtri(n+z_up,l,:),arr(n+z_up)) 
     70             CALL gradq(n,n+t_ldown,n+t_rdown,n+z_down,qi(:,l),gradtri(n+z_down,l,:),arr(n+z_down)) 
     71          END DO 
     72       END DO 
     73    END DO 
     74 
     75    !      Do l =1,llm 
     76    DO j=jj_begin,jj_end 
     77       DO i=ii_begin,ii_end 
     78          n=(j-1)*iim+i 
     79          gradq3d(n,:,:) = gradtri(n+z_up,:,:) + gradtri(n+z_down,:,:)   +   & 
     80               gradtri(n+z_rup,:,:) +  gradtri(n+z_ldown,:,:)   +   &  
     81               gradtri(n+z_lup,:,:)+    gradtri(n+z_rdown,:,:)   
     82          ar = arr(n+z_up)+arr(n+z_lup)+arr(n+z_ldown)+arr(n+z_down)+arr(n+z_rdown)+arr(n+z_rup) 
     83          gradq3d(n,:,:) = gradq3d(n,:,:)/(ar+1.e-50)  
     84       END DO 
     85    END DO 
     86    !    END DO 
     87 
     88    !============================================================================================= LIMITING  
     89    ! GO TO 120  
     90    DO l =1,llm 
     91       DO j=jj_begin,jj_end 
    11292          DO i=ii_begin,ii_end 
    11393             n=(j-1)*iim+i 
    114              gradq3d(n,:,:) = gradtri(n+z_up,:,:) + gradtri(n+z_down,:,:)   +   & 
    115                            gradtri(n+z_rup,:,:) +  gradtri(n+z_ldown,:,:)   +   &  
    116                           gradtri(n+z_lup,:,:)+    gradtri(n+z_rdown,:,:)   
    117              ar = arr(n+z_up)+arr(n+z_lup)+arr(n+z_ldown)+arr(n+z_down)+arr(n+z_rdown)+arr(n+z_rup) 
    118              gradq3d(n,:,:) = gradq3d(n,:,:)/(ar+1.e-50)  
    119            END DO 
    120          END DO     
    121 !        END DO 
    122   
    123 !============================================================================================= LIMITING  
    124 ! GO TO 120  
    125   Do l =1,llm 
    126    DO j=jj_begin,jj_end 
    127      DO i=ii_begin,ii_end 
    128        n=(j-1)*iim+i 
    129        maggrd =  dot_product(gradq3d(n,l,:),gradq3d(n,l,:)) 
    130        maggrd = sqrt(maggrd)  
    131        leng = max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), & 
    132        sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2),  & 
    133        sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2)) 
    134        maxq_c = qi(n,l) + maggrd*sqrt(leng) 
    135        minq_c = qi(n,l) - maggrd*sqrt(leng) 
    136        maxq = max(qi(n,l),qi(n+t_right,l),qi(n+t_lup,l),qi(n+t_rup,l),qi(n+t_left,l), & 
    137              qi(n+t_rdown,l),qi(n+t_ldown,l)) 
    138        minq = min(qi(n,l),qi(n+t_right,l),qi(n+t_lup,l),qi(n+t_rup,l),qi(n+t_left,l), & 
    139              qi(n+t_rdown,l),qi(n+t_ldown,l)) 
    140        alphamx = (maxq - qi(n,l)) ; alphamx = alphamx/(maxq_c - qi(n,l) ) 
    141        alphamx = max(alphamx,0.0) 
    142        alphami = (minq - qi(n,l)); alphami = alphami/(minq_c - qi(n,l)) 
    143        alphami = max(alphami,0.0)  
    144        alpha   = min(alphamx,alphami,1.0) 
    145        gradq3d(n,l,:) = alpha*gradq3d(n,l,:) 
    146       END DO 
    147     END DO 
    148    END DO 
    149          END SUBROUTINE ADVECT1           
    150 !===================================================================================================      
    151          SUBROUTINE advect2(qi,him,ue,he,bigt) 
    152   USE domain_mod 
    153   USE dimensions 
    154   USE geometry 
    155   USE metric 
    156   IMPLICIT NONE 
     94             maggrd =  dot_product(gradq3d(n,l,:),gradq3d(n,l,:)) 
     95             maggrd = sqrt(maggrd)  
     96             leng = max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), & 
     97                  sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2),  & 
     98                  sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2)) 
     99             maxq_c = qi(n,l) + maggrd*sqrt(leng) 
     100             minq_c = qi(n,l) - maggrd*sqrt(leng) 
     101             maxq = max(qi(n,l),qi(n+t_right,l),qi(n+t_lup,l),qi(n+t_rup,l),qi(n+t_left,l), & 
     102                  qi(n+t_rdown,l),qi(n+t_ldown,l)) 
     103             minq = min(qi(n,l),qi(n+t_right,l),qi(n+t_lup,l),qi(n+t_rup,l),qi(n+t_left,l), & 
     104                  qi(n+t_rdown,l),qi(n+t_ldown,l)) 
     105             alphamx = (maxq - qi(n,l)) ; alphamx = alphamx/(maxq_c - qi(n,l) ) 
     106             alphamx = max(alphamx,0.0) 
     107             alphami = (minq - qi(n,l)); alphami = alphami/(minq_c - qi(n,l)) 
     108             alphami = max(alphami,0.0)  
     109             alpha   = min(alphamx,alphami,1.0) 
     110             gradq3d(n,l,:) = alpha*gradq3d(n,l,:) 
     111          END DO 
     112       END DO 
     113    END DO 
     114  END SUBROUTINE compute_gradq3d 
     115 
     116  !===================================================================================================    
     117  SUBROUTINE compute_advect_horiz(normal,tangent,qi,gradq3d,him,ue,he,bigt) 
     118    USE domain_mod 
     119    USE dimensions 
     120    USE geometry 
     121    USE metric 
     122    IMPLICIT NONE 
     123    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3) 
     124    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3) 
    157125    REAL(rstd),INTENT(INOUT) :: qi(iim*jjm,llm) 
     126    REAL(rstd),INTENT(IN)    :: gradq3d(iim*jjm,llm,3)  
    158127    REAL(rstd),INTENT(INOUT) :: him(iim*jjm,llm) 
    159     REAL(rstd),INTENT(IN) :: ue(iim*3*jjm,llm) 
    160     REAL(rstd),INTENT(IN) :: he(3*iim*jjm,llm) 
    161     REAL(rstd),INTENT(IN)::bigt 
     128    REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm) 
     129    REAL(rstd),INTENT(IN)    :: he(3*iim*jjm,llm) ! mass flux 
     130    REAL(rstd),INTENT(IN)    :: bigt 
     131 
    162132    REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm)  
    163133    REAL(rstd) :: cc(3*iim*jjm,llm,3)  
    164134    REAL(rstd) :: v_e(3*iim*jjm,llm,3) 
    165135    REAL(rstd) :: up_e 
    166      
     136 
    167137    REAL(rstd) :: qe(3*iim*jjm,llm) 
    168138    REAL(rstd) :: ed(3)     
    … …  
    172142 
    173143 
    174 !go to  123 
    175    DO l = 1,llm 
     144    !go to  123 
     145    DO l = 1,llm 
     146       DO j=jj_begin,jj_end 
     147          DO i=ii_begin,ii_end 
     148             n=(j-1)*iim+i 
     149 
     150             up_e =1/de(n+u_right)*(                                                   & 
     151                  wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
     152                  wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
     153                  wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+                      & 
     154                  wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                    & 
     155                  wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    &  
     156                  wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+    & 
     157                  wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+    & 
     158                  wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+    & 
     159                  wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+        & 
     160                  wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         & 
     161                  ) 
     162             v_e(n+u_right,l,:)= ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:) 
     163 
     164             up_e=1/de(n+u_lup)*(                                                        & 
     165                  wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+                        & 
     166                  wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                      & 
     167                  wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                      & 
     168                  wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+                      & 
     169                  wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+                          &  
     170                  wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+          & 
     171                  wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+              & 
     172                  wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+              & 
     173                  wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+            & 
     174                  wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l)           & 
     175                  ) 
     176             v_e(n+u_lup,l,:)= ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:) 
     177 
     178             up_e=1/de(n+u_ldown)*(                                                    & 
     179                  wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    & 
     180                  wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+                    & 
     181                  wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
     182                  wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
     183                  wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+                      &  
     184                  wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+        & 
     185                  wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+      & 
     186                  wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+    & 
     187                  wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+    & 
     188                  wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l)     & 
     189                  ) 
     190             v_e(n+u_ldown,l,:)= ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:) 
     191 
     192             cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e(n+u_right,l,:)*0.5*bigt    !! redge is mid point of edge i  
     193             cc(n+u_lup,l,:)   = xyz_e(n+u_lup,:)   - v_e(n+u_lup,l,:)*0.5*bigt  
     194             cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e(n+u_ldown,l,:)*0.5*bigt  
     195          ENDDO 
     196       ENDDO 
     197    END DO 
     198    !123         continue            
     199    !========================================================================================================== 
     200    DO l = 1,llm 
     201       DO j=jj_begin-1,jj_end+1 
     202          DO i=ii_begin-1,ii_end+1 
     203             n=(j-1)*iim+i 
     204             if (ne(n,right)*ue(n+u_right,l)>0) then  
     205                ed = cc(n+u_right,l,:) - xyz_i(n,:) 
     206                qe(n+u_right,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)  
     207             else 
     208                ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:) 
     209                qe(n+u_right,l)=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed)  
     210             endif 
     211             if (ne(n,lup)*ue(n+u_lup,l)>0) then  
     212                ed = cc(n+u_lup,l,:) - xyz_i(n,:) 
     213                qe(n+u_lup,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed) 
     214             else 
     215                ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:) 
     216                qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed)  
     217             endif 
     218             if (ne(n,ldown)*ue(n+u_ldown,l)>0) then  
     219                ed = cc(n+u_ldown,l,:) - xyz_i(n,:) 
     220                qe(n+u_ldown,l)=qi(n,l)+ sum2(gradq3d(n,l,:),ed)  
     221             else 
     222                ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:) 
     223                qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed)  
     224             endif 
     225          end do 
     226       end do 
     227    END DO 
     228 
     229 
     230    DO j=jj_begin-1,jj_end+1 
     231       DO i=ii_begin-1,ii_end+1 
     232          n=(j-1)*iim+i  
     233          flxx(n+u_right,:) = he(n+u_right,:)*qe(n+u_right,:)*ne(n,right) 
     234          flxx(n+u_lup,:) = he(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup)  
     235          flxx(n+u_ldown,:) = he(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown) 
     236       ENDDO 
     237    ENDDO 
     238 
    176239    DO j=jj_begin,jj_end 
    177240       DO i=ii_begin,ii_end 
    178         n=(j-1)*iim+i 
    179  
    180                     up_e =1/de(n+u_right)*(                                                   & 
    181                          wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
    182                          wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
    183                          wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+                      & 
    184                          wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                    & 
    185                          wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    &  
    186                          wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+    & 
    187                          wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+    & 
    188                          wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+    & 
    189                          wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+        & 
    190                          wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         & 
    191                                         ) 
    192        v_e(n+u_right,l,:)= ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:) 
    193    
    194                 up_e=1/de(n+u_lup)*(                                                        & 
    195                          wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+                        & 
    196                          wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                      & 
    197                          wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                      & 
    198                          wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+                      & 
    199                          wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+                          &  
    200                          wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+          & 
    201                          wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+              & 
    202                          wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+              & 
    203                          wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+            & 
    204                          wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l)           & 
    205                                   ) 
    206        v_e(n+u_lup,l,:)= ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:) 
    207      
    208                   up_e=1/de(n+u_ldown)*(                                                    & 
    209                          wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    & 
    210                          wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+                    & 
    211                          wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
    212                          wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
    213                          wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+                      &  
    214                          wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+        & 
    215                          wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+      & 
    216                          wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+    & 
    217                          wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+    & 
    218                          wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l)     & 
    219                                       ) 
    220          v_e(n+u_ldown,l,:)= ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:) 
    221  
    222           cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e(n+u_right,l,:)*0.5*bigt    !! redge is mid point of edge i  
    223           cc(n+u_lup,l,:)   = xyz_e(n+u_lup,:)   - v_e(n+u_lup,l,:)*0.5*bigt  
    224           cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e(n+u_ldown,l,:)*0.5*bigt  
    225        ENDDO 
    226      ENDDO   
    227     END DO 
    228 !123     continue            
    229 !========================================================================================================== 
    230     DO l = 1,llm 
    231         DO j=jj_begin-1,jj_end+1 
    232          DO i=ii_begin-1,ii_end+1 
    233             n=(j-1)*iim+i 
    234         if (ne(n,right)*ue(n+u_right,l)>0) then  
    235           ed = cc(n+u_right,l,:) - xyz_i(n,:) 
    236          qe(n+u_right,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)  
    237         else 
    238            ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:) 
    239          qe(n+u_right,l)=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed)  
    240         endif 
    241         if (ne(n,lup)*ue(n+u_lup,l)>0) then  
    242            ed = cc(n+u_lup,l,:) - xyz_i(n,:) 
    243          qe(n+u_lup,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed) 
    244         else 
    245            ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:) 
    246          qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed)  
    247         endif 
    248         if (ne(n,ldown)*ue(n+u_ldown,l)>0) then  
    249           ed = cc(n+u_ldown,l,:) - xyz_i(n,:) 
    250          qe(n+u_ldown,l)=qi(n,l)+ sum2(gradq3d(n,l,:),ed)  
    251         else 
    252          ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:) 
    253         qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed)  
    254         endif 
    255         end do 
    256     end do 
    257    END DO 
    258  
    259  
    260          DO j=jj_begin-1,jj_end+1 
    261            DO i=ii_begin-1,ii_end+1 
    262             n=(j-1)*iim+i  
    263             flxx(n+u_right,:) = he(n+u_right,:)*qe(n+u_right,:)*ne(n,right) 
    264             flxx(n+u_lup,:) = he(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup)  
    265             flxx(n+u_ldown,:) = he(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown) 
    266            ENDDO 
    267          ENDDO 
    268                                  
    269      DO j=jj_begin,jj_end 
    270        DO i=ii_begin,ii_end 
    271         n=(j-1)*iim+i   
    272  
    273          dhi(n,:)= -(1/Ai(n))*(he(n+u_right,:)*ne(n,right) + he(n+u_lup,:)*ne(n,lup) & 
    274                        + he(n+u_ldown,:)*ne(n,ldown) + he(n+u_rup,:)*ne(n,rup) & 
    275                        + he(n+u_left,:)*ne(n,left) +  he(n+u_rdown,:)*ne(n,rdown) )  
    276  
    277          dqi(n,:)= -(1/Ai(n))*(flxx(n+u_right,:)+flxx(n+u_lup,:)       &  
    278                          +flxx(n+u_ldown,:) - flxx(n+u_rup,:)    & 
    279                          - flxx(n+u_left,:) - flxx(n+u_rdown,:) )    
    280          him_old(:) = him(n,:)  
    281          him(n,:) = him(n,:) + dhi(n,:)  
    282          qi(n,:) = (qi(n,:)*him_old(:) + dqi(n,:))/him(n,:)   
    283                         
    284        END DO 
    285      END DO 
    286  
    287         CONTAINS 
    288 !==================================================================================== 
    289           REAL FUNCTION sum2(a1,a2) 
    290         IMPLICIT NONE 
    291         REAL,INTENT(IN):: a1(3), a2(3) 
    292                 sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3) 
    293            END FUNCTION sum2 
    294   
    295                 END SUBROUTINE advect2  
    296 !========================================================================== 
     241          n=(j-1)*iim+i   
     242 
     243          dhi(n,:)= -(1/Ai(n))*(he(n+u_right,:)*ne(n,right) + he(n+u_lup,:)*ne(n,lup) & 
     244               + he(n+u_ldown,:)*ne(n,ldown) + he(n+u_rup,:)*ne(n,rup) & 
     245               + he(n+u_left,:)*ne(n,left) +  he(n+u_rdown,:)*ne(n,rdown) )  
     246 
     247          dqi(n,:)= -(1/Ai(n))*(flxx(n+u_right,:)+flxx(n+u_lup,:)       &  
     248               +flxx(n+u_ldown,:) - flxx(n+u_rup,:)    & 
     249               - flxx(n+u_left,:) - flxx(n+u_rdown,:) )    
     250          him_old(:) = him(n,:)  
     251          him(n,:) = him(n,:) + dhi(n,:)  
     252          qi(n,:) = (qi(n,:)*him_old(:) + dqi(n,:))/him(n,:)   
     253 
     254       END DO 
     255    END DO 
     256 
     257  CONTAINS 
     258    !==================================================================================== 
     259    REAL FUNCTION sum2(a1,a2) 
     260      IMPLICIT NONE 
     261      REAL,INTENT(IN):: a1(3), a2(3) 
     262      sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3) 
     263    END FUNCTION sum2 
     264 
     265  END SUBROUTINE compute_advect_horiz 
     266  !========================================================================== 
    297267  SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det) 
    298   USE geometry 
    299   USE metric 
    300   USE dimensions 
    301   IMPLICIT NONE 
     268    USE geometry 
     269    USE metric 
     270    USE dimensions 
     271    IMPLICIT NONE 
    302272    INTEGER, INTENT(IN) :: n0,n1,n2,n3 
    303273    REAL,INTENT(IN)     :: q(iim*jjm) 
    … …  
    313283    A(2,1)=xyz_i(n2,1) - xyz_i(n0,1); A(2,2)=xyz_i(n2,2) - xyz_i(n0,2); A(2,3)=xyz_i(n2,3) - xyz_i(n0,3)  
    314284    A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);             A(3,3)= xyz_v(n3,3) 
    315     
    316    
     285 
     286 
    317287    dq(1) = q(n1)-q(n0) 
    318288    dq(2) = q(n2)-q(n0) 
    319289    dq(3) = 0.0 
    320 !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info) 
    321         CALL determinant(A(:,1),A(:,2),A(:,3),det) 
    322         CALL determinant(dq,A(:,2),A(:,3),detx) 
    323         CALL determinant(A(:,1),dq,A(:,3),dety) 
    324         CALL determinant(A(:,1),A(:,2),dq,detz) 
    325         dq(1) = detx 
    326         dq(2) = dety 
    327         dq(3) = detz  
     290    !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info) 
     291    CALL determinant(A(:,1),A(:,2),A(:,3),det) 
     292    CALL determinant(dq,A(:,2),A(:,3),detx) 
     293    CALL determinant(A(:,1),dq,A(:,3),dety) 
     294    CALL determinant(A(:,1),A(:,2),dq,detz) 
     295    dq(1) = detx 
     296    dq(2) = dety 
     297    dq(3) = detz  
    328298  END SUBROUTINE gradq 
    329 !========================================================================== 
    330           SUBROUTINE determinant(a1,a2,a3,det) 
    331         IMPLICIT NONE 
    332         REAL, DIMENSION(3) :: a1, a2,a3 
    333         REAL ::  det,x1,x2,x3 
    334         x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2)) 
    335         x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1)) 
    336         x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1)) 
    337         det =  x1 - x2 + x3 
    338         END SUBROUTINE     
    339         END MODULE   
     299  !========================================================================== 
     300  SUBROUTINE determinant(a1,a2,a3,det) 
     301    IMPLICIT NONE 
     302    REAL, DIMENSION(3) :: a1, a2,a3 
     303    REAL ::  det,x1,x2,x3 
     304    x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2)) 
     305    x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1)) 
     306    x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1)) 
     307    det =  x1 - x2 + x3 
     308  END SUBROUTINE determinant 
     309 
     310END MODULE advect_mod 

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

    r20 r30  
    44  INTEGER,PARAMETER::iapp_tracvl= 3 
    55  REAL(rstd),SAVE :: dt 
    6    
     6 
     7  TYPE(t_field),POINTER :: f_normal(:) 
     8  TYPE(t_field),POINTER :: f_tangent(:) 
     9  TYPE(t_field),POINTER :: f_gradq3d(:) 
     10 
    711  PUBLIC init_advect_tracer, advect_tracer 
    812 
    913CONTAINS 
    10    
     14 
    1115  SUBROUTINE init_advect_tracer(dt_in) 
    12   USE advect_mod 
    13   IMPLICIT NONE 
     16    USE advect_mod 
     17    IMPLICIT NONE 
    1418    REAL(rstd),INTENT(IN) :: dt_in 
    15      
     19    REAL(rstd),POINTER :: tangent(:,:) 
     20    REAL(rstd),POINTER :: normal(:,:) 
     21    INTEGER :: ind 
     22 
    1623    dt=dt_in 
    17      
    18     CALL init_advect 
    19      
     24    CALL allocate_field(f_normal,field_u,type_real,3) 
     25    CALL allocate_field(f_tangent,field_u,type_real,3) 
     26    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3) 
     27 
     28    DO ind=1,ndomain 
     29       CALL swap_dimensions(ind) 
     30       CALL swap_geometry(ind) 
     31       normal=f_normal(ind) 
     32       tangent=f_tangent(ind) 
     33       CALL init_advect(normal,tangent) 
     34    END DO 
     35 
    2036  END SUBROUTINE init_advect_tracer 
    21    
    22    
    23   SUBROUTINE advect_tracer(f_ps,f_u, f_q) 
    24   USE icosa 
    25   USE advect_mod 
    26   USE disvert_mod 
    27   IMPLICIT NONE 
    28   TYPE(t_field),POINTER :: f_ps(:) 
    29   TYPE(t_field),POINTER :: f_u(:) 
    30   TYPE(t_field),POINTER :: f_q(:) 
    31   REAL(rstd),POINTER :: q(:,:,:) 
    32   REAL(rstd),POINTER :: u(:,:)  
    33   REAL(rstd),POINTER :: ps(:)  
    34   REAL(rstd),POINTER :: massflx(:,:) 
    35   REAL(rstd),POINTER :: rhodz(:,:) 
    36   TYPE(t_field),POINTER,SAVE :: f_massflxc(:) 
    37   TYPE(t_field),POINTER,SAVE :: f_massflx(:) 
    38   TYPE(t_field),POINTER,SAVE :: f_uc(:) 
    39   TYPE(t_field),POINTER,SAVE :: f_rhodzm1(:) 
    40   TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
    41   REAL(rstd),POINTER,SAVE :: massflxc(:,:) 
    42   REAL(rstd),POINTER,SAVE :: uc(:,:) 
    43   REAL(rstd),POINTER,SAVE :: rhodzm1(:,:) 
    44   REAL(rstd):: bigt  
    45   INTEGER :: ind,it,iapptrac,i,j,l,ij 
    46   INTEGER,SAVE :: iadvtr=0 
    47   LOGICAL,SAVE:: first=.TRUE.  
     37 
     38  SUBROUTINE advect_tracer(f_ps,f_u,f_q) 
     39    USE icosa 
     40    USE advect_mod 
     41    USE disvert_mod 
     42    IMPLICIT NONE 
     43    TYPE(t_field),POINTER :: f_ps(:) 
     44    TYPE(t_field),POINTER :: f_u(:) 
     45    TYPE(t_field),POINTER :: f_q(:) 
     46    REAL(rstd),POINTER :: q(:,:,:) 
     47    REAL(rstd),POINTER :: u(:,:)  
     48    REAL(rstd),POINTER :: ps(:)  
     49    REAL(rstd),POINTER :: massflx(:,:) 
     50    REAL(rstd),POINTER :: rhodz(:,:) 
     51    TYPE(t_field),POINTER,SAVE :: f_massflxc(:) 
     52    TYPE(t_field),POINTER,SAVE :: f_massflx(:) 
     53    TYPE(t_field),POINTER,SAVE :: f_uc(:) 
     54    TYPE(t_field),POINTER,SAVE :: f_rhodzm1(:) 
     55    TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
     56    REAL(rstd),POINTER,SAVE :: massflxc(:,:) 
     57    REAL(rstd),POINTER,SAVE :: uc(:,:) 
     58    REAL(rstd),POINTER,SAVE :: rhodzm1(:,:) 
     59    REAL(rstd):: bigt 
     60    INTEGER :: ind,it,i,j,l,ij 
     61    INTEGER,SAVE :: iadvtr=0 
     62    LOGICAL,SAVE:: first=.TRUE.  
    4863!------------------------------------------------------sarvesh 
     64       CALL transfert_request(f_ps,req_i1) 
     65       CALL transfert_request(f_u,req_e1) 
     66       CALL transfert_request(f_u,req_e1) 
     67       CALL transfert_request(f_q,req_i1) 
     68         
    4969    IF ( first ) THEN  
    50       CALL allocate_field(f_rhodz,field_t,type_real,llm)  
    51       CALL allocate_field(f_rhodzm1,field_t,type_real,llm)  
    52       CALL allocate_field(f_massflxc,field_u,type_real,llm)  
    53       CALL allocate_field(f_massflx,field_u,type_real,llm)  
    54       CALL allocate_field(f_uc,field_u,type_real,llm)  
    55       first = .FALSE. 
    56     END IF  
    57  
    58     DO ind=1,ndomain 
    59       CALL swap_dimensions(ind) 
    60       CALL swap_geometry(ind) 
    61       rhodz=f_rhodz(ind) 
    62       massflx=f_massflx(ind) 
    63       ps=f_ps(ind) 
    64       u=f_u(ind) 
    65        
    66       DO l = 1, llm 
    67         DO j=jj_begin-1,jj_end+1 
    68           DO i=ii_begin-1,ii_end+1 
    69             ij=(j-1)*iim+i 
    70             rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g  
    71           ENDDO 
    72         ENDDO 
    73       ENDDO 
    74  
    75       DO l = 1, llm 
    76         DO j=jj_begin-1,jj_end+1 
    77           DO i=ii_begin-1,ii_end+1 
    78             ij=(j-1)*iim+i 
    79             massflx(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
    80             massflx(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
    81             massflx(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
    82           ENDDO 
    83         ENDDO 
    84       ENDDO 
    85      
    86     ENDDO 
    87      
    88  
    89      
    90         IF ( iadvtr == 0 ) THEN  
    91          DO ind=1,ndomain          
    92          CALL swap_dimensions(ind) 
    93          CALL swap_geometry(ind) 
    94            rhodz=f_rhodz(ind)   
    95            rhodzm1 = f_rhodzm1(ind)  
    96            massflxc = f_massflxc(ind)  
    97            rhodzm1 = rhodz 
    98            massflxc = 0.0  
    99            uc = f_uc(ind) 
    100            uc = 0.0  
    101          END DO 
    102          CALL transfert_request(f_rhodzm1,req_i1)   !----> 
    103          CALL transfert_request(f_massflxc,req_e1) !----> 
    104          CALL transfert_request(f_massflxc,req_e1) !------> 
    105          CALL transfert_request(f_uc,req_e1) !----> 
    106          CALL transfert_request(f_uc,req_e1)  
    107         END IF 
    108  
    109         iadvtr = iadvtr + 1  
    110         iapptrac = iadvtr   
    111  
    112       DO ind=1,ndomain 
    113        CALL swap_dimensions(ind) 
    114        CALL swap_geometry(ind) 
    115             massflx=f_massflx(ind) 
    116             rhodzm1 = f_rhodzm1(ind)  
    117             massflxc = f_massflxc(ind)  
    118             massflxc = massflxc + massflx  
    119              uc = f_uc(ind)  
    120              u  = f_u(ind)  
    121              uc = uc + u  
     70       CALL allocate_field(f_rhodz,field_t,type_real,llm)  
     71       CALL allocate_field(f_rhodzm1,field_t,type_real,llm)  
     72       CALL allocate_field(f_massflxc,field_u,type_real,llm)  
     73       CALL allocate_field(f_massflx,field_u,type_real,llm)  
     74       CALL allocate_field(f_uc,field_u,type_real,llm)  
     75       first = .FALSE. 
     76    END IF 
     77 
     78    DO ind=1,ndomain 
     79       CALL swap_dimensions(ind) 
     80       CALL swap_geometry(ind) 
     81       rhodz=f_rhodz(ind) 
     82       massflx=f_massflx(ind) 
     83       ps=f_ps(ind) 
     84       u=f_u(ind) 
     85 
     86       DO l = 1, llm 
     87          DO j=jj_begin-1,jj_end+1 
     88             DO i=ii_begin-1,ii_end+1 
     89                ij=(j-1)*iim+i 
     90                rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g  
     91             ENDDO 
     92          ENDDO 
     93       ENDDO 
     94 
     95       DO l = 1, llm 
     96          DO j=jj_begin-1,jj_end+1 
     97             DO i=ii_begin-1,ii_end+1 
     98                ij=(j-1)*iim+i 
     99                massflx(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
     100                massflx(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
     101                massflx(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
     102             ENDDO 
     103          ENDDO 
     104       ENDDO 
     105    ENDDO 
     106 
     107    IF ( iadvtr == 0 ) THEN  
     108       DO ind=1,ndomain          
     109          CALL swap_dimensions(ind) 
     110          CALL swap_geometry(ind) 
     111          rhodz=f_rhodz(ind)   
     112          rhodzm1 = f_rhodzm1(ind)  
     113          massflxc = f_massflxc(ind)  ! accumulated mass fluxes 
     114          uc = f_uc(ind)              ! time-integrated normal velocity to compute back-trajectory (Miura) 
     115          rhodzm1 = rhodz 
     116          massflxc = 0.0  
     117          uc = 0.0  
    122118       END DO 
    123  
    124      IF ( iadvtr == iapp_tracvl ) THEN  
    125          bigt = dt*iapp_tracvl 
    126          DO ind=1,ndomain 
    127            CALL swap_dimensions(ind) 
    128            CALL swap_geometry(ind) 
    129               uc = f_uc(ind)  
    130               uc = uc/real(iapp_tracvl)  
    131          END DO 
     119       CALL transfert_request(f_rhodzm1,req_i1)   !----> 
     120       CALL transfert_request(f_massflxc,req_e1) !----> 
     121       CALL transfert_request(f_massflxc,req_e1) !------> 
     122       CALL transfert_request(f_uc,req_e1) !----> 
     123       CALL transfert_request(f_uc,req_e1)  
     124    END IF 
     125 
     126    iadvtr = iadvtr + 1  
     127 
     128    DO ind=1,ndomain 
     129       CALL swap_dimensions(ind) 
     130       CALL swap_geometry(ind) 
     131       massflx  = f_massflx(ind) 
     132       massflxc = f_massflxc(ind)  
     133       uc = f_uc(ind)  
     134       u  = f_u(ind)  
     135       massflxc = massflxc + massflx  
     136       uc = uc + u  
     137    END DO 
     138 
     139    IF ( iadvtr == iapp_tracvl ) THEN  
     140       PRINT *, 'Advection scheme' 
     141       bigt = dt*iapp_tracvl 
     142       DO ind=1,ndomain 
     143          CALL swap_dimensions(ind) 
     144          CALL swap_geometry(ind) 
     145          uc = f_uc(ind)  
     146          uc = uc/real(iapp_tracvl)  
     147          massflxc = f_massflxc(ind) 
     148          massflxc = massflxc*dt 
     149          ! now massflx is time-integrated 
     150       END DO 
    132151 
    133152       CALL vlsplt(f_q,f_rhodzm1,f_massflxc,2.0,f_uc,bigt)  
    134             iadvtr = 0 
    135      END IF 
     153       iadvtr = 0 
     154    END IF 
    136155  END SUBROUTINE advect_tracer 
    137    
    138 !==============================================================================   
    139   SUBROUTINE advtrac(massflx,wgg) 
    140   USE domain_mod 
    141   USE dimensions 
    142   USE grid_param 
    143   USE geometry 
    144   USE metric 
    145   USE disvert_mod 
    146   IMPLICIT NONE 
    147     REAL(rstd),INTENT(IN) :: massflx(iim*3*jjm,llm) 
    148     REAL(rstd),INTENT(OUT) :: wgg(iim*jjm,llm) 
    149       
    150     INTEGER :: i,j,ij,l 
    151     REAL(rstd) :: convm(iim*jjm,llm)  
    152     
    153      DO l = 1, llm 
    154       DO j=jj_begin,jj_end 
    155        DO i=ii_begin,ii_end 
    156         ij=(j-1)*iim+i 
    157         convm(ij,l)= 1/(Ai(ij))*(ne(ij,right)*massflx(ij+u_right,l)   +  & 
    158                                 ne(ij,rup)*massflx(ij+u_rup,l)       +  & 
    159                                 ne(ij,lup)*massflx(ij+u_lup,l)       +  & 
    160                                 ne(ij,left)*massflx(ij+u_left,l)     +  & 
    161                                 ne(ij,ldown)*massflx(ij+u_ldown,l)   +  & 
    162                                 ne(ij,rdown)*massflx(ij+u_rdown,l)) 
    163         ENDDO 
    164       ENDDO 
    165      ENDDO 
    166  
    167     DO  l = llm-1, 1, -1 
    168      DO j=jj_begin,jj_end 
    169       DO i=ii_begin,ii_end 
    170         ij=(j-1)*iim+i 
    171         convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
    172       ENDDO 
    173      ENDDO 
    174     ENDDO 
    175  
    176     !!! Compute vertical velocity 
    177   DO l = 1,llm-1 
    178     DO j=jj_begin,jj_end 
    179       DO i=ii_begin,ii_end 
    180         ij=(j-1)*iim+i 
    181         wgg( ij, l+1 ) = (convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ))  
    182       ENDDO 
    183     ENDDO 
    184   ENDDO 
    185  
    186     DO j=jj_begin,jj_end 
    187      DO i=ii_begin,ii_end 
    188       ij=(j-1)*iim+i 
    189       wgg(ij,1)  = 0. 
    190      ENDDO 
    191    ENDDO 
    192  END SUBROUTINE advtrac  
    193  
    194         SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt)  
    195   USE field_mod 
    196   USE domain_mod 
    197   USE dimensions 
    198   USE grid_param 
    199   USE geometry 
    200   USE metric 
    201   USE advect_mod 
    202   IMPLICIT NONE 
    203  
    204   TYPE(t_field),POINTER :: f_q(:) 
    205   TYPE(t_field),POINTER :: f_u(:) 
    206   TYPE(t_field),POINTER :: f_rhodz(:)  
    207   TYPE(t_field),POINTER :: f_massflx(:) 
    208   TYPE(t_field),POINTER,SAVE :: f_wg(:) 
    209   TYPE(t_field),POINTER,SAVE :: f_zm(:) 
    210   TYPE(t_field),POINTER,SAVE :: f_zq(:) 
    211  
    212   REAL(rstd)::bigt,dt 
    213   REAL(rstd),POINTER :: q(:,:,:) 
    214   REAL(rstd),POINTER :: u(:,:) 
    215   REAL(rstd),POINTER :: rhodz(:,:)  
    216   REAL(rstd),POINTER :: massflx(:,:)  
    217   REAL(rstd),POINTER,SAVE :: wg(:,:) 
    218   REAL(rstd),POINTER,SAVE::zq(:,:,:)  
    219   REAL(rstd),POINTER,SAVE::zm(:,:)  
    220  
    221   REAL(rstd)::pente_max 
    222   LOGICAL,SAVE::first = .true.  
    223   INTEGER :: i,ij,l,j,ind,k 
    224   REAL(rstd) :: zzpbar, zzw  
    225   REAL::qvmax,qvmin  
    226    
    227      IF ( first ) THEN  
     156 
     157  SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt) 
     158    USE field_mod 
     159    USE domain_mod 
     160    USE dimensions 
     161    USE grid_param 
     162    USE geometry 
     163    USE metric 
     164    USE advect_mod 
     165    IMPLICIT NONE 
     166 
     167    TYPE(t_field),POINTER :: f_q(:) 
     168    TYPE(t_field),POINTER :: f_u(:) 
     169    TYPE(t_field),POINTER :: f_rhodz(:)  
     170    TYPE(t_field),POINTER :: f_massflx(:) 
     171 
     172    TYPE(t_field),POINTER,SAVE :: f_wg(:) 
     173    TYPE(t_field),POINTER,SAVE :: f_zm(:) 
     174    TYPE(t_field),POINTER,SAVE :: f_zq(:) 
     175 
     176    REAL(rstd)::bigt,dt 
     177    REAL(rstd),POINTER :: q(:,:,:) 
     178    REAL(rstd),POINTER :: u(:,:) 
     179    REAL(rstd),POINTER :: rhodz(:,:)  
     180    REAL(rstd),POINTER :: massflx(:,:)  
     181    REAL(rstd),POINTER :: wg(:,:) 
     182    REAL(rstd),POINTER :: zq(:,:,:)  
     183    REAL(rstd),POINTER :: zm(:,:)  
     184    REAL(rstd),POINTER :: tangent(:,:) 
     185    REAL(rstd),POINTER :: normal(:,:) 
     186    REAL(rstd),POINTER :: gradq3d(:,:,:) 
     187 
     188    REAL(rstd)::pente_max 
     189    LOGICAL,SAVE::first = .true.  
     190    INTEGER :: i,ij,l,j,ind,k 
     191    REAL(rstd) :: zzpbar, zzw  
     192    REAL::qvmax,qvmin  
     193 
     194    IF ( first ) THEN  
    228195       CALL allocate_field(f_wg,field_t,type_real,llm) 
    229196       CALL allocate_field(f_zm,field_t,type_real,llm) 
    230197       CALL allocate_field(f_zq,field_t,type_real,llm,nqtot)  
    231198       first = .FALSE. 
    232       END IF  
    233  
    234      DO ind=1,ndomain 
    235        CALL swap_dimensions(ind) 
    236        CALL swap_geometry(ind) 
    237          q=f_q(ind)  
    238          rhodz=f_rhodz(ind)  
    239          zq=f_zq(ind) 
    240          zm=f_zm(ind)  
    241          zm = rhodz ; zq = q   
    242          wg = f_wg(ind) 
    243          wg = 0.0 
    244          massflx=f_massflx(ind)  
    245        CALL advtrac(massflx,wg)  
    246      END DO 
    247  
    248 !   CALL transfert_request(f_wg,req_i1)  
    249  
    250       DO ind=1,ndomain 
    251         CALL swap_dimensions(ind) 
    252         CALL swap_geometry(ind) 
    253          zq=f_zq(ind) 
    254          zm=f_zm(ind)  
    255          wg=f_wg(ind) 
    256          wg=wg*0.5 
     199    END IF 
     200 
     201    DO ind=1,ndomain 
     202       CALL swap_dimensions(ind) 
     203       CALL swap_geometry(ind) 
     204       q=f_q(ind)  
     205       rhodz=f_rhodz(ind)  
     206       zq=f_zq(ind) 
     207       zm=f_zm(ind)  
     208       zm = rhodz ; zq = q   
     209       wg = f_wg(ind) 
     210       wg = 0.0 
     211       massflx=f_massflx(ind)  
     212       CALL advtrac(massflx,wg) ! compute vertical mass fluxes 
     213    END DO 
     214 
     215    DO ind=1,ndomain 
     216       CALL swap_dimensions(ind) 
     217       CALL swap_geometry(ind) 
     218       zq=f_zq(ind) 
     219       zm=f_zm(ind)  
     220       wg=f_wg(ind) 
     221       wg=wg*0.5 
    257222       DO k = 1, nqtot 
    258         CALL vlz(zq(:,:,k),2.,zm,wg) 
     223          CALL vlz(zq(:,:,k),2.,zm,wg) 
    259224       END DO 
    260       END DO 
    261  
    262         DO ind=1,ndomain 
    263           CALL swap_dimensions(ind) 
    264           CALL swap_geometry(ind) 
    265           CALL swap_advect(ind) 
    266           zq=f_zq(ind) 
    267           zq = f_zq(ind) 
    268           zm = f_zm(ind) 
    269           massflx =f_massflx(ind) 
    270           u = f_u(ind)  
    271          DO k = 1,nqtot 
    272            CALL advect1(zq(:,:,k)) 
    273            CALL advect2(zq(:,:,k),zm,u,massflx,bigt)  
    274          END DO 
    275         END DO 
    276  
    277         DO ind=1,ndomain 
    278            CALL swap_dimensions(ind) 
    279            CALL swap_geometry(ind) 
    280             q = f_q(ind) 
    281            zq = f_zq(ind)  
    282            zm = f_zm(ind)  
    283            wg = f_wg(ind) 
    284           DO k = 1,nqtot 
    285             CALL vlz(zq(:,:,k),2.,zm,wg) 
    286           END DO  
    287             q = zq  
    288         END DO 
    289  
    290   END SUBROUTINE vlsplt  
    291  
    292          
    293      SUBROUTINE vlz(q,pente_max,masse,wgg) 
    294 !c 
    295 !c     Auteurs:   P.Le Van, F.Hourdin, F.Forget  
    296 !c 
    297 !c    ******************************************************************** 
    298 !c     Shema  d'advection " pseudo amont " . 
    299 !c    ******************************************************************** 
    300       USE icosa 
    301       IMPLICIT NONE 
    302 !c 
    303 !c   Arguments: 
    304 !c   ---------- 
    305       REAL masse(iim*jjm,llm),pente_max 
    306       REAL q(iim*jjm,llm) 
    307       REAL wgg(iim*jjm,llm),w(iim*jjm,llm+1)  
    308       REAL dq(iim*jjm,llm) 
    309       INTEGER i,ij,l,j,ii 
    310 !c 
    311       REAL wq(iim*jjm,llm+1),newmasse 
    312  
    313       REAL dzq(iim*jjm,llm),dzqw(iim*jjm,llm),adzqw(iim*jjm,llm),dzqmax 
    314       REAL sigw 
    315  
    316       REAL      SSUM 
    317  
    318  
    319          w(:,1:llm) = wgg(:,:) 
    320          w(:,llm+1) = 0.0  
    321  
    322 !c    On oriente tout dans le sens de la pression c'est a dire dans le 
    323 !c    sens de W 
    324  
    325        DO l=2,llm 
    326         DO j=jj_begin,jj_end 
    327          DO i=ii_begin,ii_end 
    328             ij=(j-1)*iim+i 
    329             dzqw(ij,l)=q(ij,l-1)-q(ij,l) 
    330             adzqw(ij,l)=abs(dzqw(ij,l)) 
    331          ENDDO 
    332         ENDDO 
    333        ENDDO 
    334  
    335        DO l=2,llm-1 
    336         DO j=jj_begin,jj_end 
    337          DO i=ii_begin,ii_end 
    338             ij=(j-1)*iim+i 
    339            IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 
     225    END DO 
     226 
     227    DO ind=1,ndomain 
     228       CALL swap_dimensions(ind) 
     229       CALL swap_geometry(ind) 
     230       zq=f_zq(ind) 
     231       zq = f_zq(ind) 
     232       zm = f_zm(ind) 
     233       massflx =f_massflx(ind) 
     234       u = f_u(ind) 
     235       tangent = f_tangent(ind) 
     236       normal = f_normal(ind) 
     237       gradq3d = f_gradq3d(ind) 
     238 
     239       DO k = 1,nqtot 
     240          CALL compute_gradq3d(zq(:,:,k),gradq3d) 
     241          CALL compute_advect_horiz(tangent,normal,zq(:,:,k),gradq3d,zm,u,massflx,bigt)  
     242       END DO 
     243    END DO 
     244 
     245    DO ind=1,ndomain 
     246       CALL swap_dimensions(ind) 
     247       CALL swap_geometry(ind) 
     248       q = f_q(ind) 
     249       zq = f_zq(ind)  
     250       zm = f_zm(ind)  
     251       wg = f_wg(ind) 
     252       DO k = 1,nqtot 
     253          CALL vlz(zq(:,:,k),2.,zm,wg) 
     254       END DO 
     255       q = zq  
     256    END DO 
     257 
     258  END SUBROUTINE vlsplt 
     259 
     260  !==============================================================================   
     261  SUBROUTINE advtrac(massflx,wgg) 
     262    USE domain_mod 
     263    USE dimensions 
     264    USE grid_param 
     265    USE geometry 
     266    USE metric 
     267    USE disvert_mod 
     268    IMPLICIT NONE 
     269    REAL(rstd),INTENT(IN) :: massflx(iim*3*jjm,llm) 
     270    REAL(rstd),INTENT(OUT) :: wgg(iim*jjm,llm) 
     271 
     272    INTEGER :: i,j,ij,l 
     273    REAL(rstd) :: convm(iim*jjm,llm)  
     274 
     275    DO l = 1, llm 
     276       DO j=jj_begin,jj_end 
     277          DO i=ii_begin,ii_end 
     278             ij=(j-1)*iim+i 
     279             ! divergence of horizontal flux 
     280             convm(ij,l)= 1/(Ai(ij))*(ne(ij,right)*massflx(ij+u_right,l)   +  & 
     281                  ne(ij,rup)*massflx(ij+u_rup,l)       +  & 
     282                  ne(ij,lup)*massflx(ij+u_lup,l)       +  & 
     283                  ne(ij,left)*massflx(ij+u_left,l)     +  & 
     284                  ne(ij,ldown)*massflx(ij+u_ldown,l)   +  & 
     285                  ne(ij,rdown)*massflx(ij+u_rdown,l)) 
     286          ENDDO 
     287       ENDDO 
     288    ENDDO 
     289 
     290    ! accumulate divergence from top of model 
     291    DO  l = llm-1, 1, -1 
     292       DO j=jj_begin,jj_end 
     293          DO i=ii_begin,ii_end 
     294             ij=(j-1)*iim+i 
     295             convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
     296          ENDDO 
     297       ENDDO 
     298    ENDDO 
     299 
     300!!! Compute vertical velocity 
     301    DO l = 1,llm-1 
     302       DO j=jj_begin,jj_end 
     303          DO i=ii_begin,ii_end 
     304             ij=(j-1)*iim+i 
     305             wgg( ij, l+1 ) = (convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ))  
     306          ENDDO 
     307       ENDDO 
     308    ENDDO 
     309 
     310    DO j=jj_begin,jj_end 
     311       DO i=ii_begin,ii_end 
     312          ij=(j-1)*iim+i 
     313          wgg(ij,1)  = 0. 
     314       ENDDO 
     315    ENDDO 
     316  END SUBROUTINE advtrac 
     317 
     318  SUBROUTINE vlz(q,pente_max,masse,wgg) 
     319    !c 
     320    !c     Auteurs:   P.Le Van, F.Hourdin, F.Forget  
     321    !c 
     322    !c    ******************************************************************** 
     323    !c     Shema  d'advection " pseudo amont " . 
     324    !c    ******************************************************************** 
     325    USE icosa 
     326    IMPLICIT NONE 
     327    !c 
     328    !c   Arguments: 
     329    !c   ---------- 
     330    REAL masse(iim*jjm,llm),pente_max 
     331    REAL q(iim*jjm,llm) 
     332    REAL wgg(iim*jjm,llm),w(iim*jjm,llm+1)  
     333    REAL dq(iim*jjm,llm) 
     334    INTEGER i,ij,l,j,ii 
     335    !c 
     336    REAL wq(iim*jjm,llm+1),newmasse 
     337 
     338    REAL dzq(iim*jjm,llm),dzqw(iim*jjm,llm),adzqw(iim*jjm,llm),dzqmax 
     339    REAL sigw 
     340 
     341    REAL      SSUM 
     342 
     343 
     344    w(:,1:llm) = -wgg(:,:)  ! w>0 for downward transport 
     345    w(:,llm+1) = 0.0  
     346 
     347    !c    On oriente tout dans le sens de la pression c'est a dire dans le 
     348    !c    sens de W 
     349 
     350    DO l=2,llm 
     351       DO j=jj_begin,jj_end 
     352          DO i=ii_begin,ii_end 
     353             ij=(j-1)*iim+i 
     354             dzqw(ij,l)=q(ij,l-1)-q(ij,l) 
     355             adzqw(ij,l)=abs(dzqw(ij,l)) 
     356          ENDDO 
     357       ENDDO 
     358    ENDDO 
     359 
     360    DO l=2,llm-1 
     361       DO j=jj_begin,jj_end 
     362          DO i=ii_begin,ii_end 
     363             ij=(j-1)*iim+i 
     364             IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 
    340365                dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1)) 
    341            ELSE 
     366             ELSE 
    342367                dzq(ij,l)=0. 
    343            ENDIF 
    344             dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1)) 
    345             dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l)) 
    346          ENDDO 
    347         ENDDO 
    348         ENDDO 
    349  
    350        DO l=2,llm-1 
    351         DO j=jj_begin,jj_end 
    352          DO i=ii_begin,ii_end 
    353             ij=(j-1)*iim+i 
    354             dzq(ij,1)=0. 
    355            dzq(ij,llm)=0. 
    356          ENDDO 
    357         ENDDO 
    358       ENDDO 
    359 !c --------------------------------------------------------------- 
    360 !c   .... calcul des termes d'advection verticale  ....... 
    361 !c --------------------------------------------------------------- 
    362  
    363 !c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq 
    364  
    365         DO l = 1,llm-1 
    366            DO j=jj_begin,jj_end 
    367            DO i=ii_begin,ii_end 
     368             ENDIF 
     369             dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1)) 
     370             dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l)) 
     371          ENDDO 
     372       ENDDO 
     373    ENDDO 
     374 
     375    DO l=2,llm-1 
     376       DO j=jj_begin,jj_end 
     377          DO i=ii_begin,ii_end 
     378             ij=(j-1)*iim+i 
     379             dzq(ij,1)=0. 
     380             dzq(ij,llm)=0. 
     381          ENDDO 
     382      ENDDO 
     383    ENDDO 
     384    !c --------------------------------------------------------------- 
     385    !c   .... calcul des termes d'advection verticale  ....... 
     386    !c --------------------------------------------------------------- 
     387 
     388    !c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq 
     389 
     390    DO l = 1,llm-1 
     391       DO j=jj_begin,jj_end 
     392          DO i=ii_begin,ii_end 
    368393             ij=(j-1)*iim+i 
    369394             IF(w(ij,l+1).gt.0.) THEN 
    370              sigw=w(ij,l+1)/masse(ij,l+1) 
    371              wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1)) 
     395                ! upwind only if downward transport 
     396                sigw=w(ij,l+1)/masse(ij,l+1) 
     397                wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1)) 
    372398             ELSE 
    373              sigw=w(ij,l+1)/masse(ij,l) 
    374              wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l)) 
    375             ENDIF 
    376            ENDDO 
    377          ENDDO 
    378          END DO 
    379  
    380           DO j=jj_begin,jj_end 
    381            DO i=ii_begin,ii_end 
    382              ij=(j-1)*iim+i 
    383              wq(ij,llm+1)=0. 
    384              wq(ij,1)=0. 
    385            ENDDO 
    386           END DO 
    387  
    388       DO l=1,llm 
    389        DO j=jj_begin,jj_end 
    390          DO i=ii_begin,ii_end 
    391              ij=(j-1)*iim+i 
    392             newmasse=masse(ij,l)+(w(ij,l+1)-w(ij,l)) 
    393             dq(ij,l) = (wq(ij,l+1)-wq(ij,l)) !================>>>> 
    394             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))/newmasse     
    395                           
    396             masse(ij,l)=newmasse 
    397            ENDDO 
    398         ENDDO 
    399       END DO     
    400       RETURN 
    401       END SUBROUTINE vlz 
     399                ! upwind only if upward transport 
     400                sigw=w(ij,l+1)/masse(ij,l) 
     401                wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l)) 
     402             ENDIF 
     403          ENDDO 
     404       ENDDO 
     405    END DO 
     406 
     407    DO j=jj_begin,jj_end 
     408       DO i=ii_begin,ii_end 
     409          ij=(j-1)*iim+i 
     410          wq(ij,llm+1)=0. 
     411          wq(ij,1)=0. 
     412       ENDDO 
     413    END DO 
     414 
     415    DO l=1,llm 
     416       DO j=jj_begin,jj_end 
     417          DO i=ii_begin,ii_end 
     418             ij=(j-1)*iim+i 
     419             ! masse -= dw/dz but w>0 <=> downward 
     420             newmasse=masse(ij,l)+(w(ij,l+1)-w(ij,l))  
     421!             dq(ij,l) = (wq(ij,l+1)-wq(ij,l)) !================>>>> 
     422             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))/newmasse     
     423!             masse(ij,l)=newmasse 
     424          ENDDO 
     425       ENDDO 
     426    END DO 
     427    RETURN 
     428  END SUBROUTINE vlz 
    402429 
    403430END MODULE advect_tracer_mod 

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

    r20 r30  
    112112  REAL(rstd),POINTER :: dtheta_rhodz(:,:) 
    113113  REAL(rstd),POINTER :: du(:,:) 
    114   INTEGER :: ind 
     114  INTEGER :: ind,ij 
    115115 
    116116   
    … …  
    119119    CALL transfert_request(f_theta_rhodz,req_i1)  
    120120    CALL transfert_request(f_u,req_e1) 
    121     CALL transfert_request(f_u,req_e1)  
     121!    CALL transfert_request(f_u,req_e1)  
    122122    
    123123 
    … …  
    136136      dtheta_rhodz=f_dtheta_rhodz(ind) 
    137137      du=f_du(ind) 
    138        
     138!      ij=(jj_end-1-1)*iim+ii_begin 
     139!      PRINT *,"--> ind=",ind,ij 
     140!      PRINT *,u(ij+u_right,1) 
     141!      PRINT *,u(ij+u_rup,1) 
     142!      PRINT *,u(ij+u_lup,1) 
     143!      PRINT *,u(ij+u_left,1) 
     144!      PRINT *,u(ij+u_ldown,1) 
     145!      PRINT *,u(ij+u_rdown,1) 
     146 
     147!      ij=(jj_end-1-1)*iim+ii_end 
     148!      PRINT *,"--> ind=",ind,ij 
     149!      PRINT *,u(ij+u_right,1) 
     150!      PRINT *,u(ij+u_rup,1) 
     151!      PRINT *,u(ij+u_lup,1) 
     152!      PRINT *,u(ij+u_left,1) 
     153!      PRINT *,u(ij+u_ldown,1) 
     154!      PRINT *,u(ij+u_rdown,1)       
    139155!$OMP PARALLEL DEFAULT(SHARED) 
    140156      CALL compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du) 
    … …  
    143159 
    144160    CALL transfert_request(f_out_u,req_e1) 
    145     CALL transfert_request(f_out_u,req_e1)  
     161!    CALL transfert_request(f_out_u,req_e1)  
    146162 
    147163!    CALL vorticity(f_u,f_out_z) 
    148 !    CALL kinetic(f_du,f_out) 
    149164 
    150165    IF (mod(it,itau_out)==0 ) THEN 
    151166      CALL writefield("ps",f_ps) 
    152 !      CALL writefield("dps",f_dps) 
     167      CALL writefield("dps",f_dps) 
    153168!      CALL writefield("theta_rhodz",f_theta_rhodz) 
     169!    CALL kinetic(f_u,f_out) 
     170!      CALL writefield("Ki",f_out) 
    154171!      CALL writefield("dtheta_rhodz",f_dtheta_rhodz) 
    155172      CALL vorticity(f_u,f_out_z) 
    … …  
    163180!        CALL writefield("Ki",f_out,ind) 
    164181!        CALL writefield("vort",f_out_z,ind) 
     182!        CALL writefield("dps",f_dps,ind) 
    165183!      ENDDO 
    166184 
    … …  
    383401    
    384402!!!  Compute mass flux 
    385 !! question à thomas : meilleure pondération de la masse sur les liens ? 
     403!! question ï¿œ thomas : meilleure pondï¿œration de la masse sur les liens ? 
    386404 
    387405  DO l = 1, llm 
    … …  
    420438        ij=(j-1)*iim+i   
    421439! signe ? attention d (rho theta dz) 
     440        ! dtheta_rhodz =  -div(flux.theta) 
    422441        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  & 
    423442                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  &   
    … …  
    440459      DO i=ii_begin,ii_end 
    441460        ij=(j-1)*iim+i 
    442 !signe ? 
     461        ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
    443462        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  & 
    444463                                ne(ij,rup)*Fe(ij+u_rup,l)       +  &   
    … …  
    480499    DO i=ii_begin,ii_end 
    481500      ij=(j-1)*iim+i 
     501      ! dps/dt = -int(div flux)dz 
    482502      dps(ij)=-convm(ij,1) * g  
    483503    ENDDO 
    … …  
    492512      DO i=ii_begin,ii_end 
    493513        ij=(j-1)*iim+i 
     514        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 
     515        ! => w>0 for upward transport 
    494516        w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ) 
    495517      ENDDO 
    … …  
    664686    DO j=jj_begin,jj_end 
    665687      DO i=ii_begin,ii_end 
     688         ! ww>0 <=> upward transport 
    666689        ij=(j-1)*iim+i 
    667690        ww            =   0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) ) 
    668         dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww  
     691        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww 
    669692        dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww 
    670693      ENDDO 

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

    r20 r30  
    2424  TYPE(t_request),POINTER :: req(:)  
    2525  TYPE(t_request),POINTER :: req_u(:)    
    26   PUBLIC :: allocate_caldyn,swap_caldyn,init_williamson,caldyn,write_caldyn,Compute_PV,Compute_enstrophy 
     26  PUBLIC :: allocate_caldyn,swap_caldyn,caldyn,write_caldyn,Compute_PV,Compute_enstrophy 
    2727CONTAINS 
    2828 

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

    r20 r30  
    22  USE icosa 
    33 
    4   TYPE(t_field),POINTER,SAVE :: f_gradrot(:) 
    5   TYPE(t_request),POINTER :: req_dissip(:)  
    6   TYPE(t_field),POINTER,SAVE :: f_due(:) 
     4  PRIVATE 
     5 
     6  TYPE(t_field),POINTER,SAVE :: f_due_diss1(:) 
     7  TYPE(t_field),POINTER,SAVE :: f_due_diss2(:) 
     8 
     9  TYPE(t_field),POINTER,SAVE :: f_theta(:) 
     10  TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:) 
     11  TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:) 
     12   
    713  
    8   INTEGER,PARAMETER :: nitergdiv=1 
    9   INTEGER,PARAMETER :: nitergrot=1 
    10   INTEGER,PARAMETER :: niterdivgrad=1 
     14  INTEGER,SAVE :: nitergdiv=1 
     15  INTEGER,SAVE :: nitergrot=1 
     16  INTEGER,SAVE :: niterdivgrad=1 
    1117 
    1218  REAL,ALLOCATABLE,SAVE :: tau_graddiv(:) 
    … …  
    1824  REAL,SAVE :: cdivgrad 
    1925   
    20   INTEGER :: idissip 
    21   REAL    :: dtdissip 
    22    
     26  INTEGER,SAVE :: idissip 
     27  REAL,SAVE    :: dtdissip 
     28   
     29  PUBLIC init_dissip, dissip 
    2330CONTAINS 
    2431 
    … …  
    2633  USE icosa 
    2734  IMPLICIT NONE   
    28     CALL allocate_field(f_gradrot,field_u,type_real,llm) 
    29     CALL allocate_field(f_due,field_u,type_real,llm) 
     35    CALL allocate_field(f_due_diss1,field_u,type_real,llm) 
     36    CALL allocate_field(f_due_diss2,field_u,type_real,llm) 
     37    CALL allocate_field(f_theta,field_t,type_real,llm) 
     38    CALL allocate_field(f_dtheta_diss,field_t,type_real,llm) 
     39    CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm) 
     40     
    3041    ALLOCATE(tau_graddiv(llm)) 
    3142    ALLOCATE(tau_gradrot(llm))     
    … …  
    3647  USE icosa 
    3748  USE disvert_mod 
     49  USE mpi_mod 
     50  USE mpipara 
    3851   
    3952  IMPLICIT NONE 
    … …  
    5669   
    5770             
    58   INTEGER :: i,j,n,ind,it 
     71  INTEGER :: i,j,n,ind,it,iter 
    5972 
    6073    CALL allocate_dissip 
    … …  
    6780    CALL getin("tau_graddiv",tau) 
    6881    tau_graddiv(:)=tau 
     82 
     83    CALL getin("nitergdiv",nitergdiv) 
    6984     
    7085    tau_gradrot(:)=5000 
    71     CALL getin("tau_gradrot",tau) 
     86    CALL getin("tau_gradrot",tau_gradrot) 
    7287    tau_gradrot(:)=tau 
     88 
     89    CALL getin("nitergrot",nitergrot) 
    7390     
    7491 
    … …  
    7693    CALL getin("tau_divgrad",tau) 
    7794    tau_divgrad(:)=tau 
    78    
    79     CALL create_request(field_u,req_dissip) 
    80  
    81     DO ind=1,ndomain 
    82       DO i=ii_begin,ii_end 
    83         CALL request_add_point(ind,i,jj_begin-1,req_dissip,rup) 
    84         CALL request_add_point(ind,i+1,jj_begin-1,req_dissip,lup) 
    85       ENDDO 
    86  
    87       DO j=jj_begin,jj_end 
    88         CALL request_add_point(ind,ii_end+1,j,req_dissip,left) 
    89         CALL request_add_point(ind,ii_end+1,j-1,req_dissip,lup) 
    90       ENDDO     
    91      
    92       DO i=ii_begin,ii_end 
    93         CALL request_add_point(ind,i,jj_end+1,req_dissip,ldown) 
    94         CALL request_add_point(ind,i-1,jj_end+1,req_dissip,rdown) 
    95       ENDDO     
    96  
    97       DO j=jj_begin,jj_end 
    98         CALL request_add_point(ind,ii_begin-1,j,req_dissip,right) 
    99         CALL request_add_point(ind,ii_begin-1,j+1,req_dissip,rdown) 
    100       ENDDO    
    101  
    102       DO i=ii_begin+1,ii_end-1 
    103         CALL request_add_point(ind,i,jj_begin,req_dissip,right) 
    104         CALL request_add_point(ind,i,jj_end,req_dissip,right) 
    105       ENDDO 
    106      
    107       DO j=jj_begin+1,jj_end-1 
    108         CALL request_add_point(ind,ii_begin,j,req_dissip,rup) 
    109         CALL request_add_point(ind,ii_end,j,req_dissip,rup) 
    110       ENDDO    
    111  
    112       CALL request_add_point(ind,ii_begin+1,jj_begin,req_dissip,left) 
    113       CALL request_add_point(ind,ii_begin,jj_begin+1,req_dissip,ldown) 
    114       CALL request_add_point(ind,ii_begin+1,jj_end,req_dissip,left) 
    115       CALL request_add_point(ind,ii_end,jj_begin+1,req_dissip,ldown) 
    116       
    117     ENDDO 
     95 
     96    CALL getin("niterdivgrad",niterdivgrad) 
     97   
     98!   CALL create_request(field_u,req_dissip) 
     99 
     100!   DO ind=1,ndomain 
     101!     DO i=ii_begin,ii_end 
     102!       CALL request_add_point(ind,i,jj_begin-1,req_dissip,rup) 
     103!       CALL request_add_point(ind,i+1,jj_begin-1,req_dissip,lup) 
     104!     ENDDO 
     105 
     106!     DO j=jj_begin,jj_end 
     107!       CALL request_add_point(ind,ii_end+1,j,req_dissip,left) 
     108!       CALL request_add_point(ind,ii_end+1,j-1,req_dissip,lup) 
     109!     ENDDO     
     110!    
     111!     DO i=ii_begin,ii_end 
     112!       CALL request_add_point(ind,i,jj_end+1,req_dissip,ldown) 
     113!       CALL request_add_point(ind,i-1,jj_end+1,req_dissip,rdown) 
     114!     ENDDO     
     115 
     116!     DO j=jj_begin,jj_end 
     117!       CALL request_add_point(ind,ii_begin-1,j,req_dissip,right) 
     118!       CALL request_add_point(ind,ii_begin-1,j+1,req_dissip,rdown) 
     119!     ENDDO    
     120 
     121!     DO i=ii_begin+1,ii_end-1 
     122!       CALL request_add_point(ind,i,jj_begin,req_dissip,right) 
     123!       CALL request_add_point(ind,i,jj_end,req_dissip,right) 
     124!     ENDDO 
     125!    
     126!     DO j=jj_begin+1,jj_end-1 
     127!       CALL request_add_point(ind,ii_begin,j,req_dissip,rup) 
     128!       CALL request_add_point(ind,ii_end,j,req_dissip,rup) 
     129!     ENDDO    
     130 
     131!     CALL request_add_point(ind,ii_begin+1,jj_begin,req_dissip,left) 
     132!     CALL request_add_point(ind,ii_begin,jj_begin+1,req_dissip,ldown) 
     133!     CALL request_add_point(ind,ii_begin+1,jj_end,req_dissip,left) 
     134!     CALL request_add_point(ind,ii_end,jj_begin+1,req_dissip,ldown) 
     135!     
     136!   ENDDO 
    118137 
    119138 
    … …  
    145164 
    146165 
    147     DO it=1,500 
    148  
    149       CALL transfert_request(f_u,req_dissip) 
    150       CALL transfert_request(f_u,req_dissip) 
    151  
     166    DO it=1,20 
     167       
    152168      dumax=0 
    153       DO ind=1,ndomain 
    154         CALL swap_dimensions(ind) 
    155         CALL swap_geometry(ind) 
    156         u=f_u(ind) 
    157         du=f_du(ind) 
    158         CALL gradiv(u,du,1) 
    159       ENDDO 
    160       CALL transfert_request(f_du,req_dissip) 
    161       CALL transfert_request(f_du,req_dissip) 
     169      DO iter=1,nitergdiv 
     170        CALL transfert_request(f_u,req_e1) 
     171        DO ind=1,ndomain 
     172          CALL swap_dimensions(ind) 
     173          CALL swap_geometry(ind) 
     174          u=f_u(ind) 
     175          du=f_du(ind) 
     176          CALL compute_gradiv(u,du,1) 
     177          u=du 
     178        ENDDO 
     179      ENDDO 
     180 
     181      CALL transfert_request(f_du,req_e1) 
    162182       
    163183      DO ind=1,ndomain 
    … …  
    166186        u=f_u(ind) 
    167187        du=f_du(ind) 
    168          
     188          
    169189        DO j=jj_begin,jj_end 
    170190          DO i=ii_begin,ii_end 
    … …  
    177197      ENDDO 
    178198 
     199      IF (using_mpi) CALL MPI_ALLREDUCE(dumax,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 
     200 
    179201      DO ind=1,ndomain 
    180202        CALL swap_dimensions(ind) 
    … …  
    212234 
    213235 
    214     DO it=1,500 
     236    DO it=1,20 
    215237  
    216       CALL transfert_request(f_u,req_dissip) 
    217       CALL transfert_request(f_u,req_dissip) 
    218  
    219238      dumax=0 
    220       DO ind=1,ndomain 
    221         CALL swap_dimensions(ind) 
    222         CALL swap_geometry(ind) 
    223         u=f_u(ind) 
    224         du=f_du(ind) 
    225         CALL gradrot(u,du,1) 
    226       ENDDO 
    227  
    228       CALL transfert_request(f_du,req_dissip) 
    229       CALL transfert_request(f_du,req_dissip) 
     239      DO iter=1,nitergrot 
     240        CALL transfert_request(f_u,req_e1) 
     241        DO ind=1,ndomain 
     242          CALL swap_dimensions(ind) 
     243          CALL swap_geometry(ind) 
     244          u=f_u(ind) 
     245          du=f_du(ind) 
     246          CALL compute_gradrot(u,du,1) 
     247          u=du 
     248        ENDDO 
     249      ENDDO 
     250 
     251      CALL transfert_request(f_du,req_e1) 
    230252       
    231253      DO ind=1,ndomain 
    … …  
    244266        ENDDO 
    245267      ENDDO 
     268 
     269      IF (using_mpi) CALL MPI_ALLREDUCE(dumax,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 
    246270       
    247271      DO ind=1,ndomain 
    … …  
    253277      ENDDO 
    254278       
    255 !      PRINT *,"gradrot : it :",it ,": dumax",(dumax+dumaxm1)/2 
    256279      PRINT *,"gradrot : it :",it ,": dumax",dumax 
    257280 
    … …  
    259282    PRINT *,"gradrot : dumax",dumax 
    260283   
    261     cgradrot=dumax**(-1/nitergrot)  
     284    cgradrot=dumax**(-1./nitergrot)  
    262285    PRINT *,"cgradrot : ",cgradrot 
    263286    
    … …  
    279302    ENDDO 
    280303 
    281     DO it=1,500 
    282  
    283       CALL transfert_request(f_theta,req_i1) 
     304    DO it=1,20 
    284305 
    285306      dthetamax=0 
    286       DO ind=1,ndomain 
    287         CALL swap_dimensions(ind) 
    288         CALL swap_geometry(ind) 
    289         theta=f_theta(ind) 
    290         dtheta=f_dtheta(ind) 
    291         CALL divgrad(theta,dtheta,1) 
     307      DO iter=1,niterdivgrad 
     308        CALL transfert_request(f_theta,req_i1) 
     309        DO ind=1,ndomain 
     310          CALL swap_dimensions(ind) 
     311          CALL swap_geometry(ind) 
     312          theta=f_theta(ind) 
     313          dtheta=f_dtheta(ind) 
     314          CALL compute_divgrad(theta,dtheta,1) 
     315          theta=dtheta 
     316        ENDDO 
    292317      ENDDO 
    293318!      CALL writefield("divgrad",f_dtheta) 
    … …  
    308333        ENDDO 
    309334      ENDDO 
    310  
     335      IF (using_mpi) CALL MPI_ALLREDUCE(dthetamax,dthetamax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 
    311336      PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax 
    312337 
    … …  
    323348    PRINT *,"divgrad : divgrad",dthetamax 
    324349   
    325     cdivgrad=dthetamax**(-1/niterdivgrad)  
     350    cdivgrad=dthetamax**(-1./niterdivgrad)  
    326351    PRINT *,"cdivgrad : ",cdivgrad 
    327352     
    … …  
    358383    TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    359384    TYPE(t_field),POINTER :: f_dtheta_rhodz(:) 
    360     REAL(rstd),POINTER         :: ue(:,:) 
     385 
    361386    REAL(rstd),POINTER         :: due(:,:) 
    362     REAL(rstd),POINTER         :: ps(:) 
    363     REAL(rstd),POINTER         :: theta_rhodz(:,:) 
     387    REAL(rstd),POINTER         :: due_diss1(:,:) 
     388    REAL(rstd),POINTER         :: due_diss2(:,:) 
    364389    REAL(rstd),POINTER         :: dtheta_rhodz(:,:) 
    365  
    366     REAL(rstd)                 :: theta(iim*jjm,llm) 
    367     REAL(rstd)                 :: du_dissip(3*iim*jjm,llm) 
    368     REAL(rstd)                 :: dtheta_dissip(iim*jjm,llm) 
    369     REAL(rstd)                 :: dtheta_rhodz_dissip(iim*jjm,llm) 
     390    REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:) 
    370391 
    371392    INTEGER :: ind 
    372393    INTEGER :: l,i,j,n 
    373394     
    374     CALL transfert_request(f_ue,req_e1) 
    375     CALL transfert_request(f_ue,req_e1) 
    376     CALL transfert_request(f_theta_rhodz,req_i1) 
    377     CALL transfert_request(f_ps,req_i1) 
    378  
     395    CALL gradiv(f_ue,f_due_diss1) 
     396    CALL gradrot(f_ue,f_due_diss2) 
     397    CALL theta_rhodz2theta(f_ps,f_theta_rhodz,f_theta) 
     398    CALL divgrad(f_theta,f_dtheta_diss) 
     399    CALL theta2theta_rhodz(f_ps,f_dtheta_diss,f_dtheta_rhodz_diss) 
     400     
     401    DO ind=1,ndomain 
     402      CALL swap_dimensions(ind) 
     403      CALL swap_geometry(ind) 
     404      due=f_due(ind)  
     405      due_diss1=f_due_diss1(ind) 
     406      due_diss2=f_due_diss2(ind) 
     407      dtheta_rhodz=f_dtheta_rhodz(ind) 
     408      dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind) 
     409             
     410      DO l=1,llm 
     411        DO j=jj_begin,jj_end 
     412          DO i=ii_begin,ii_end 
     413            n=(j-1)*iim+i 
     414 
     415            due(n+u_right,l)=due(n+u_right,l)-0.5*(due_diss1(n+u_right,l)/tau_graddiv(l)+due_diss2(n+u_right,l)/tau_gradrot(l)) 
     416            due(n+u_lup,l)=due(n+u_lup,l)    -0.5*(due_diss1(n+u_lup,l)  /tau_graddiv(l)+due_diss2(n+u_lup,l)  /tau_gradrot(l)) 
     417            due(n+u_ldown,l)=due(n+u_ldown,l)-0.5*(due_diss1(n+u_ldown,l)/tau_graddiv(l)+due_diss2(n+u_ldown,l)/tau_gradrot(l)) 
     418 
     419            dtheta_rhodz(n,l)=dtheta_rhodz(n,l)-0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l) 
     420 
     421          ENDDO 
     422        ENDDO 
     423      ENDDO 
     424    ENDDO 
     425       
     426  END SUBROUTINE dissip 
     427 
     428  SUBROUTINE gradiv(f_ue,f_due) 
     429  USE icosa 
     430  IMPLICIT NONE 
     431    TYPE(t_field),POINTER :: f_ue(:) 
     432    TYPE(t_field),POINTER :: f_due(:) 
     433    REAL(rstd),POINTER    :: ue(:,:) 
     434    REAL(rstd),POINTER    :: due(:,:) 
     435    INTEGER :: ind 
     436    INTEGER :: it 
     437        
    379438    DO ind=1,ndomain 
    380439      CALL swap_dimensions(ind) 
    … …  
    382441      ue=f_ue(ind) 
    383442      due=f_due(ind)  
    384       ps=f_ps(ind)  
    385       theta_rhodz=f_theta_rhodz(ind) 
    386       dtheta_rhodz=f_dtheta_rhodz(ind)  
    387  
    388 !$OMP PARALLEL DEFAULT(SHARED) 
    389      CALL compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz) 
    390 !$OMP END PARALLEL 
    391        
    392     ENDDO 
    393        
    394   END SUBROUTINE dissip 
    395  
    396  
    397   SUBROUTINE compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz) 
    398   USE icosa 
    399   USE theta2theta_rhodz_mod 
     443      due=ue 
     444    ENDDO 
     445 
     446    DO it=1,nitergdiv 
     447         
     448      CALL transfert_request(f_due,req_e1) 
     449         
     450      DO ind=1,ndomain 
     451        CALL swap_dimensions(ind) 
     452        CALL swap_geometry(ind) 
     453        due=f_due(ind)  
     454        CALL compute_gradiv(due,due,llm) 
     455      ENDDO 
     456    ENDDO 
     457 
     458  END SUBROUTINE gradiv 
     459   
     460 
     461  SUBROUTINE gradrot(f_ue,f_due) 
     462  USE icosa 
    400463  IMPLICIT NONE 
    401     REAL(rstd)         :: ue(3*iim*jjm,llm) 
    402     REAL(rstd)         :: due(3*iim*jjm,llm) 
    403     REAL(rstd)         :: ps(iim*jjm) 
    404     REAL(rstd)         :: theta_rhodz(iim*jjm,llm) 
    405     REAL(rstd)         :: dtheta_rhodz(iim*jjm,llm) 
    406  
    407     REAL(rstd),SAVE,ALLOCATABLE :: theta(:,:) 
    408     REAL(rstd),SAVE,ALLOCATABLE :: du_dissip(:,:) 
    409     REAL(rstd),SAVE,ALLOCATABLE :: dtheta_dissip(:,:) 
    410     REAL(rstd),SAVE,ALLOCATABLE :: dtheta_rhodz_dissip(:,:) 
    411  
     464    TYPE(t_field),POINTER :: f_ue(:) 
     465    TYPE(t_field),POINTER :: f_due(:) 
     466    REAL(rstd),POINTER    :: ue(:,:) 
     467    REAL(rstd),POINTER    :: due(:,:) 
    412468    INTEGER :: ind 
    413     INTEGER :: l,i,j,n 
    414  
    415 !$OMP BARRIER 
    416 !$OMP MASTER 
    417       ALLOCATE(theta(iim*jjm,llm)) 
    418       ALLOCATE(du_dissip(3*iim*jjm,llm)) 
    419       ALLOCATE(dtheta_dissip(iim*jjm,llm)) 
    420       ALLOCATE(dtheta_rhodz_dissip(iim*jjm,llm)) 
    421 !$OMP END MASTER 
    422 !$OMP BARRIER 
    423      
    424       CALL gradiv(ue,du_dissip,llm) 
    425       DO l=1,llm 
    426 !$OMP DO 
    427         DO j=jj_begin,jj_end 
    428           DO i=ii_begin,ii_end 
    429             n=(j-1)*iim+i 
    430             due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_graddiv(l)*0.5 
    431             due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_graddiv(l)*0.5 
    432             due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_graddiv(l)*0.5 
    433           ENDDO 
    434         ENDDO 
    435       ENDDO 
    436        
    437       CALL gradrot(ue,du_dissip,llm) 
    438  
    439       DO l=1,llm 
    440 !$OMP DO 
    441         DO j=jj_begin,jj_end 
    442           DO i=ii_begin,ii_end 
    443             n=(j-1)*iim+i 
    444             due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_gradrot(l)*0.5 
    445             due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_gradrot(l)*0.5 
    446             due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_gradrot(l)*0.5 
    447           ENDDO 
    448         ENDDO 
    449       ENDDO  
    450        
    451       CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1) 
    452       CALL divgrad(theta,dtheta_dissip,llm) 
    453       CALL compute_theta2theta_rhodz(ps,dtheta_dissip,dtheta_rhodz_dissip,0) 
    454  
    455       DO l=1,llm 
    456 !$OMP DO 
    457         DO j=jj_begin,jj_end 
    458           DO i=ii_begin,ii_end 
    459             n=(j-1)*iim+i 
    460             dtheta_rhodz(n,l)=dtheta_rhodz(n,l)+dtheta_rhodz_dissip(n,l)/tau_divgrad(l)*0.5 
    461           ENDDO 
    462         ENDDO 
    463       ENDDO  
    464  
    465 !$OMP BARRIER 
    466 !$OMP MASTER 
    467       DEALLOCATE(theta) 
    468       DEALLOCATE(du_dissip) 
    469       DEALLOCATE(dtheta_dissip) 
    470       DEALLOCATE(dtheta_rhodz_dissip) 
    471 !$OMP END MASTER 
    472 !$OMP BARRIER       
    473  
    474   END SUBROUTINE compute_dissip 
    475    
    476    
    477   SUBROUTINE gradiv(ue,gradivu_e,ll) 
     469    INTEGER :: it 
     470        
     471    DO ind=1,ndomain 
     472      CALL swap_dimensions(ind) 
     473      CALL swap_geometry(ind) 
     474      ue=f_ue(ind) 
     475      due=f_due(ind)  
     476      due=ue 
     477    ENDDO 
     478 
     479    DO it=1,nitergrot 
     480         
     481      CALL transfert_request(f_due,req_e1) 
     482         
     483      DO ind=1,ndomain 
     484        CALL swap_dimensions(ind) 
     485        CALL swap_geometry(ind) 
     486        due=f_due(ind)  
     487        CALL compute_gradrot(due,due,llm) 
     488      ENDDO 
     489 
     490    ENDDO 
     491 
     492  END SUBROUTINE gradrot 
     493   
     494  SUBROUTINE divgrad(f_theta,f_dtheta) 
     495  USE icosa 
     496  IMPLICIT NONE 
     497    TYPE(t_field),POINTER :: f_theta(:) 
     498    TYPE(t_field),POINTER :: f_dtheta(:) 
     499    REAL(rstd),POINTER    :: theta(:,:) 
     500    REAL(rstd),POINTER    :: dtheta(:,:) 
     501    INTEGER :: ind 
     502    INTEGER :: it 
     503        
     504    DO ind=1,ndomain 
     505      CALL swap_dimensions(ind) 
     506      CALL swap_geometry(ind) 
     507      theta=f_theta(ind) 
     508      dtheta=f_dtheta(ind)  
     509      dtheta=theta 
     510    ENDDO 
     511 
     512    DO it=1,niterdivgrad 
     513         
     514      CALL transfert_request(f_dtheta,req_i1) 
     515         
     516      DO ind=1,ndomain 
     517        CALL swap_dimensions(ind) 
     518        CALL swap_geometry(ind) 
     519        dtheta=f_dtheta(ind)  
     520        CALL compute_divgrad(dtheta,dtheta,llm) 
     521      ENDDO 
     522 
     523    ENDDO 
     524 
     525  END SUBROUTINE divgrad 
     526    
     527   
     528 
     529!  SUBROUTINE compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz) 
     530!  USE icosa 
     531!  USE theta2theta_rhodz_mod 
     532!  IMPLICIT NONE 
     533!    REAL(rstd)         :: ue(3*iim*jjm,llm) 
     534!    REAL(rstd)         :: due(3*iim*jjm,llm) 
     535!    REAL(rstd)         :: ps(iim*jjm) 
     536!    REAL(rstd)         :: theta_rhodz(iim*jjm,llm) 
     537!    REAL(rstd)         :: dtheta_rhodz(iim*jjm,llm) 
     538! 
     539!    REAL(rstd),SAVE,ALLOCATABLE :: theta(:,:) 
     540!    REAL(rstd),SAVE,ALLOCATABLE :: du_dissip(:,:) 
     541!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_dissip(:,:) 
     542!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_rhodz_dissip(:,:) 
     543! 
     544!    INTEGER :: ind 
     545!    INTEGER :: l,i,j,n 
     546! 
     547!!$OMP BARRIER 
     548!!$OMP MASTER 
     549!      ALLOCATE(theta(iim*jjm,llm)) 
     550!      ALLOCATE(du_dissip(3*iim*jjm,llm)) 
     551!      ALLOCATE(dtheta_dissip(iim*jjm,llm)) 
     552!      ALLOCATE(dtheta_rhodz_dissip(iim*jjm,llm)) 
     553!!$OMP END MASTER 
     554!!$OMP BARRIER 
     555!     
     556!      CALL gradiv(ue,du_dissip,llm) 
     557!      DO l=1,llm 
     558!!$OMP DO 
     559!        DO j=jj_begin,jj_end 
     560!          DO i=ii_begin,ii_end 
     561!            n=(j-1)*iim+i 
     562!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_graddiv(l)*0.5 
     563!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_graddiv(l)*0.5 
     564!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_graddiv(l)*0.5 
     565!          ENDDO 
     566!        ENDDO 
     567!      ENDDO 
     568!       
     569!      CALL gradrot(ue,du_dissip,llm) 
     570! 
     571!      DO l=1,llm 
     572!!$OMP DO 
     573!        DO j=jj_begin,jj_end 
     574!          DO i=ii_begin,ii_end 
     575!            n=(j-1)*iim+i 
     576!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_gradrot(l)*0.5 
     577!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_gradrot(l)*0.5 
     578!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_gradrot(l)*0.5 
     579!          ENDDO 
     580!        ENDDO 
     581!      ENDDO  
     582!       
     583!      CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1) 
     584!      CALL divgrad(theta,dtheta_dissip,llm) 
     585!      CALL compute_theta2theta_rhodz(ps,dtheta_dissip,dtheta_rhodz_dissip,0) 
     586! 
     587!      DO l=1,llm 
     588!!$OMP DO 
     589!        DO j=jj_begin,jj_end 
     590!          DO i=ii_begin,ii_end 
     591!            n=(j-1)*iim+i 
     592!            dtheta_rhodz(n,l)=dtheta_rhodz(n,l)+dtheta_rhodz_dissip(n,l)/tau_divgrad(l)*0.5 
     593!          ENDDO 
     594!        ENDDO 
     595!      ENDDO  
     596! 
     597!!$OMP BARRIER 
     598!!$OMP MASTER 
     599!      DEALLOCATE(theta) 
     600!      DEALLOCATE(du_dissip) 
     601!      DEALLOCATE(dtheta_dissip) 
     602!      DEALLOCATE(dtheta_rhodz_dissip) 
     603!!$OMP END MASTER 
     604!!$OMP BARRIER       
     605! 
     606!  END SUBROUTINE compute_dissip 
     607   
     608   
     609  SUBROUTINE compute_gradiv(ue,gradivu_e,ll) 
    478610  USE icosa 
    479611  IMPLICIT NONE 
    … …  
    528660        DO i=ii_begin,ii_end 
    529661          n=(j-1)*iim+i 
    530           gradivu_e(n+u_right,l)=gradivu_e(n+u_right,l)*cgraddiv 
    531           gradivu_e(n+u_lup,l)=gradivu_e(n+u_lup,l)*cgraddiv 
    532           gradivu_e(n+u_ldown,l)=gradivu_e(n+u_ldown,l)*cgraddiv 
     662          gradivu_e(n+u_right,l)=-gradivu_e(n+u_right,l)*cgraddiv 
     663          gradivu_e(n+u_lup,l)=-gradivu_e(n+u_lup,l)*cgraddiv 
     664          gradivu_e(n+u_ldown,l)=-gradivu_e(n+u_ldown,l)*cgraddiv 
    533665        ENDDO 
    534666      ENDDO 
    … …  
    542674 
    543675     
    544   END SUBROUTINE gradiv 
    545    
    546   SUBROUTINE divgrad(theta,divgrad_i,ll) 
     676  END SUBROUTINE compute_gradiv 
     677   
     678  SUBROUTINE compute_divgrad(theta,divgrad_i,ll) 
    547679  USE icosa 
    548680  IMPLICIT NONE 
    … …  
    598730        DO i=ii_begin,ii_end 
    599731          n=(j-1)*iim+i 
    600           divgrad_i(n,l)=divgrad_i(n,l)*cdivgrad 
     732          divgrad_i(n,l)=-divgrad_i(n,l)*cdivgrad 
    601733        ENDDO 
    602734      ENDDO 
    … …  
    609741!$OMP BARRIER 
    610742 
    611   END SUBROUTINE divgrad 
    612  
    613      
    614   SUBROUTINE gradrot(ue,gradrot_e,ll) 
     743  END SUBROUTINE compute_divgrad 
     744 
     745     
     746  SUBROUTINE compute_gradrot(ue,gradrot_e,ll) 
    615747  USE icosa 
    616748  IMPLICIT NONE 
    … …  
    668800          n=(j-1)*iim+i 
    669801     
    670           gradrot_e(n+u_right,l)=gradrot_e(n+u_right,l)*cgradrot        
    671           gradrot_e(n+u_lup,l)=gradrot_e(n+u_lup,l)*cgradrot 
    672           gradrot_e(n+u_ldown,l)=gradrot_e(n+u_ldown,l)*cgradrot 
     802          gradrot_e(n+u_right,l)=-gradrot_e(n+u_right,l)*cgradrot        
     803          gradrot_e(n+u_lup,l)=-gradrot_e(n+u_lup,l)*cgradrot 
     804          gradrot_e(n+u_ldown,l)=-gradrot_e(n+u_ldown,l)*cgradrot 
    673805         
    674806        ENDDO 
    … …  
    682814!$OMP BARRIER  
    683815    
    684   END SUBROUTINE gradrot 
     816  END SUBROUTINE compute_gradrot 
    685817 
    686818 

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

    r20 r30  
    66  REAL(rstd), SAVE, ALLOCATABLE,TARGET :: presnivs(:) 
    77 
    8 CONTAINS  
     8CONTAINS 
    99!========================================================================= 
    1010 
    … …  
    2121  END SUBROUTINE init_disvert   
    2222 
    23  
    2423  SUBROUTINE disvert(ap,bp,presnivs) 
    2524  USE icosa 
    … …  
    2827  REAL(rstd),INTENT(OUT) :: bp(:) 
    2928  REAL(rstd),INTENT(OUT) :: presnivs(:) 
     29  ! reads from run.def : ncar_dz, ncar_T0, ncar_p0, ncar_disvert_c 
     30  INTEGER :: l,cindx 
     31  REAL(rstd) :: H, eta_top, eta 
     32 
     33  ncar_dz=400 ; CALL getin('ncar_dz',ncar_dz); 
     34 
     35! SELECT CASE(  
     36! pressure profile depends on test case 
     37! coded here for 1-x (transport) 
     38 
     39  ncar_T0=300; CALL getin('ncar_T0',ncar_T0) 
     40  ncar_p0=1e5; CALL getin('ncar_p0',ncar_p0) 
     41  cindx=1 ; CALL getin('ncar_disvert_c',cindx) 
     42 
     43  H = ncar_T0*cpp*kappa/g ! height scale R.T0/g with R=cpp*kappa 
     44  eta_top = exp(-llm*ncar_dz/H) 
     45  
     46  do l = 1,llm+1 
     47     eta = exp(-(l-1)*ncar_dz/H) 
     48     bp(l) = ((eta - eta_top)/(1 - eta_top))**cindx 
     49     ap(l) = ncar_p0 * ( eta - bp(l) ) 
     50  ENDDO 
     51  bp(1)=1. 
     52  ap(1)=0. 
     53  bp(llm+1) = 0 
    3054   
    31   REAL(rstd) :: sig(llm+1) 
    32   REAL(rstd) :: sigtop 
    33   REAL(rstd),PARAMETER:: p0=100000.0  
    34   INTEGER :: l,cindx 
    35   REAL(rstd) :: hdz, ehdz  
     55  DO l = 1, llm 
     56     presnivs(l) = 0.5 *( ap(l)+bp(l)*ncar_p0 + ap(l+1)+bp(l+1)*ncar_p0 ) 
     57  ENDDO 
    3658 
     59  PRINT *, 'Vertical placement of model levels according to DCMIP Appendix E.3' 
     60  PRINT *, 'Parameters : ncar_dz=', ncar_dz, '  ncar_p0=',ncar_p0, '  ncar_disvert_c=',cindx 
     61  PRINT *, 'Isothermal amtosphere with ncar_T0=',ncar_T0  
    3762 
    38          hdz = 400.*g/(300.*287.)  
    39          ehdz = exp(-hdz)  
    40  
    41     do l = 1,llm+1 
    42         sig(l) = ehdz**(l-1)  
    43     end do   
    44          
    45        sigtop = sig(llm+1)  
    46         cindx = 1  
    47  
    48     bp(llm+1) =   0. 
    49     DO l = 1, llm 
    50       bp(l) = (sig(l) - sigtop)/(1 - sigtop) 
    51       bp(l) = bp(l)**cindx 
    52       ap(l) = p0 * ( sig(l) - bp(l) ) 
    53     ENDDO 
    54     bp(1)=1. 
    55     ap(1)=0. 
    56     ap(llm+1) = p0 * ( sig(llm+1) - bp(llm+1) ) 
    57    
    58     DO l = 1, llm 
    59       presnivs(l) = 0.5 *( ap(l)+bp(l)*p0 + ap(l+1)+bp(l+1)*p0 ) 
    60    ENDDO 
    61  
    62   END SUBROUTINE disvert 
     63END SUBROUTINE disvert 
    6364 
    6465END  MODULE disvert_ncar_mod 

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

    r20 r30  
    1919    INTEGER,POINTER  :: assign_i(:,:) 
    2020    INTEGER,POINTER  :: assign_j(:,:) 
     21    INTEGER,POINTER  :: edge_assign_domain(:,:,:) 
     22    INTEGER,POINTER  :: edge_assign_i(:,:,:) 
     23    INTEGER,POINTER  :: edge_assign_j(:,:,:) 
     24    INTEGER,POINTER  :: edge_assign_pos(:,:,:) 
    2125    REAL,POINTER     :: xyz(:,:,:) 
    2226    REAL,POINTER     :: neighbour(:,:,:,:) 
    … …  
    5559  TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain(:) 
    5660  INTEGER :: ndomain 
    57  
    58  
     61  TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain_glo(:) 
     62  INTEGER :: ndomain_glo 
     63 
     64  INTEGER,ALLOCATABLE,SAVE :: domglo_rank(:) 
     65  INTEGER,ALLOCATABLE,SAVE :: domglo_loc_ind(:) 
     66  INTEGER,ALLOCATABLE,SAVE :: domloc_glo_ind(:) 
     67   
    5968CONTAINS 
    6069 
    … …  
    6675  TYPE(t_domain),POINTER :: d 
    6776  
    68     ndomain=nsplit_i*nsplit_j*nb_face 
    69     ALLOCATE(domain(ndomain)) 
    70    
     77    ndomain_glo=nsplit_i*nsplit_j*nb_face 
     78    ALLOCATE(domain_glo(ndomain_glo)) 
     79    ALLOCATE(domglo_rank(ndomain_glo)) 
     80    ALLOCATE(domglo_loc_ind(ndomain_glo)) 
     81 
    7182    ind=0 
    7283    DO nf=1,nb_face 
    73       DO nj=1,nsplit_i 
    74         DO ni=1,nsplit_j 
     84      DO nj=1,nsplit_j 
     85        DO ni=1,nsplit_i 
    7586          ind=ind+1 
    76           d=>domain(ind) 
     87          d=>domain_glo(ind) 
    7788          quotient=(iim_glo/nsplit_i) 
    7889          rest=MOD(iim_glo,nsplit_i) 
    … …  
    8596          ENDIF 
    8697          d%ii_end_glo=d%ii_begin_glo+d%ii_nb-1 
     98  
     99          IF (ni/=nsplit_i) THEN  
     100            d%ii_nb=d%ii_nb+1 
     101            d%ii_end_glo=d%ii_end_glo+1 
     102          ENDIF 
     103          
    87104         
    88105          quotient=(jjm_glo/nsplit_j) 
    … …  
    95112            d%jj_begin_glo=(quotient+1)*rest+(nj-1-rest) * quotient+1 
    96113          ENDIF 
    97  
    98114          d%jj_end_glo=d%jj_begin_glo+d%jj_nb-1 
     115 
     116          IF (nj/=nsplit_j) THEN  
     117            d%jj_nb=d%jj_nb+1 
     118            d%jj_end_glo=d%jj_end_glo+1 
     119          ENDIF 
     120 
     121 
    99122          d%iim=d%ii_nb+2*halo 
    100123          d%jjm=d%jj_nb+2*halo 
    … …  
    107130          ALLOCATE(d%assign_i(d%iim,d%jjm)) 
    108131          ALLOCATE(d%assign_j(d%iim,d%jjm)) 
     132          ALLOCATE(d%edge_assign_domain(0:5,d%iim,d%jjm)) 
     133          ALLOCATE(d%edge_assign_i(0:5,d%iim,d%jjm)) 
     134          ALLOCATE(d%edge_assign_j(0:5,d%iim,d%jjm)) 
     135          ALLOCATE(d%edge_assign_pos(0:5,d%iim,d%jjm)) 
    109136          ALLOCATE(d%delta(d%iim,d%jjm)) 
    110137          ALLOCATE(d%xyz(3,d%iim,d%jjm)) 
    … …  
    119146  END SUBROUTINE create_domain 
    120147   
     148  SUBROUTINE copy_domain(d1,d2) 
     149  IMPLICIT NONE 
     150  INTEGER :: face 
     151  TYPE(t_domain),TARGET,INTENT(IN) :: d1 
     152  TYPE(t_domain), INTENT(OUT) :: d2 
     153   
     154    d2%iim = d1%iim 
     155    d2%jjm = d1%jjm 
     156    d2%ii_begin = d1%ii_begin 
     157    d2%jj_begin = d1%jj_begin 
     158    d2%ii_end = d1%ii_end 
     159    d2%jj_end = d1%jj_end 
     160    d2%ii_nb =  d1%ii_nb 
     161    d2%jj_nb = d1%jj_nb 
     162    d2%ii_begin_glo = d1%ii_begin_glo 
     163    d2%jj_begin_glo = d1%jj_begin_glo 
     164    d2%ii_end_glo = d1%ii_end_glo 
     165    d2%jj_end_glo = d1%jj_end_glo 
     166    d2%assign_domain => d1%assign_domain 
     167    d2%assign_i => d1%assign_i 
     168    d2%assign_j => d1%assign_j 
     169    d2%edge_assign_domain => d1%edge_assign_domain 
     170    d2%edge_assign_i => d1%edge_assign_i 
     171    d2%edge_assign_j => d1%edge_assign_j 
     172    d2%edge_assign_pos => d1%edge_assign_pos 
     173    d2%xyz => d1%xyz 
     174    d2%neighbour => d1%neighbour 
     175    d2%delta => d1%delta 
     176    d2%vertex => d1%vertex 
     177    d2%own => d1%own 
     178    d2%ne => d1%ne 
     179     
     180    d2%t_right = d1%t_right 
     181    d2%t_rup = d1%t_rup 
     182    d2%t_lup = d1%t_lup 
     183    d2%t_left = d1%t_left 
     184    d2%t_ldown = d1%t_ldown 
     185    d2%t_rdown = d1%t_rdown 
     186 
     187    d2%u_right = d1%u_right 
     188    d2%u_rup = d1%u_rup 
     189    d2%u_lup = d1%u_lup 
     190    d2%u_left = d1%u_left 
     191    d2%u_ldown = d1%u_ldown 
     192    d2%u_rdown = d1%u_rdown 
     193 
     194    d2%z_rup = d1%z_rup 
     195    d2%z_up = d1%z_up 
     196    d2%z_lup = d1%z_lup 
     197    d2%z_ldown = d1%z_ldown 
     198    d2%z_down = d1%z_down 
     199    d2%z_rdown = d1%z_rdown 
     200    
     201    d2%t_pos = d1%t_pos 
     202    d2%u_pos = d1%u_pos 
     203    d2%z_pos = d1%z_pos 
     204     
     205  END SUBROUTINE copy_domain 
     206   
    121207   
    122208  SUBROUTINE assign_cell 
    123209  USE metric 
    124210  IMPLICIT NONE 
    125     INTEGER :: ind_d,ind,ind2 
     211    INTEGER :: ind_d,ind,ind2,e 
    126212    INTEGER :: nf 
    127213    INTEGER :: i,j,k,ii,jj 
    128214    TYPE(t_domain),POINTER :: d 
     215    INTEGER :: delta 
    129216      
    130217     
    131     DO ind_d=1,ndomain 
    132       d=>domain(ind_d) 
     218    DO ind_d=1,ndomain_glo 
     219      d=>domain_glo(ind_d) 
    133220      nf=d%face 
    134221      DO j=d%jj_begin,d%jj_end 
    … …  
    137224          jj=d%jj_begin_glo-d%jj_begin+j 
    138225          ind=vertex_glo(ii,jj,nf)%ind 
     226          delta=vertex_glo(ii,jj,nf)%delta 
    139227          IF (cell_glo(ind)%assign_face==nf) THEN  
    140228            cell_glo(ind)%assign_domain=ind_d 
    … …  
    142230            cell_glo(ind)%assign_j=j 
    143231          ENDIF 
     232          IF ( i+1 <= d%ii_end ) CALL assign_edge(ind_d,ind,i,j,delta,0) 
     233          IF ( j+1 <= d%jj_end ) CALL assign_edge(ind_d,ind,i,j,delta,1) 
     234          IF ( i-1 >= d%ii_begin .AND. j+1<=d%jj_end ) CALL assign_edge(ind_d,ind,i,j,delta,2) 
     235          IF ( i-1 >= d%ii_begin ) CALL assign_edge(ind_d,ind,i,j,delta,3) 
     236          IF ( j-1 >= d%jj_begin ) CALL assign_edge(ind_d,ind,i,j,delta,4) 
     237          IF ( i+1 <= d%ii_end .AND. j-1 >=d%jj_begin ) CALL assign_edge(ind_d,ind,i,j,delta,5) 
    144238        ENDDO 
    145239      ENDDO 
    146240    ENDDO 
    147241     
    148     DO ind_d=1,ndomain 
    149       d=>domain(ind_d) 
     242     
     243    DO ind_d=1,ndomain_glo 
     244      d=>domain_glo(ind_d) 
    150245      nf=d%face 
    151246      DO j=d%jj_begin-1,d%jj_end+1 
    … …  
    157252          d%assign_i(i,j)=cell_glo(ind)%assign_i 
    158253          d%assign_j(i,j)=cell_glo(ind)%assign_j 
     254          delta=vertex_glo(ii,jj,nf)%delta 
    159255          d%delta(i,j)=vertex_glo(ii,jj,nf)%delta 
    160256          DO k=0,5 
    … …  
    162258            d%neighbour(:,k,i,j)=cell_glo(ind2)%xyz(:) 
    163259            d%ne(k,i,j)=vertex_glo(ii,jj,nf)%ne(k) 
     260            e=cell_glo(ind)%edge(MOD(k+delta+6,6)) 
     261            d%edge_assign_domain(k,i,j)=edge_glo(e)%assign_domain 
     262            d%edge_assign_i(k,i,j)=edge_glo(e)%assign_i 
     263            d%edge_assign_j(k,i,j)=edge_glo(e)%assign_j 
     264            d%edge_assign_pos(k,i,j)=edge_glo(e)%assign_pos 
    164265          ENDDO 
    165266          d%xyz(:,i,j)=cell_glo(ind)%xyz(:) 
    … …  
    171272        ENDDO 
    172273      ENDDO 
    173     ENDDO     
     274    ENDDO 
     275 
     276  CONTAINS 
     277 
     278    SUBROUTINE assign_edge(ind_d,ind,i,j,delta,k) 
     279      INTEGER :: ind_d,ind,i,j,delta,k 
     280      INTEGER :: e 
     281      e=cell_glo(ind)%edge(MOD(k+delta+6,6)) 
     282      edge_glo(e)%assign_domain=ind_d 
     283      edge_glo(e)%assign_i=i 
     284      edge_glo(e)%assign_j=j 
     285      edge_glo(e)%assign_pos=k 
     286!      PRINT*,"Assign_edge",ind_d,ind,i,j,k,e 
     287     END  SUBROUTINE assign_edge 
     288          
    174289  END SUBROUTINE assign_cell 
    175290 
    … …  
    182297    REAL(rstd) :: ng1(3),ng2(3)  
    183298 
    184     DO ind_d=1,ndomain 
    185       d=>domain(ind_d) 
     299    DO ind_d=1,ndomain_glo 
     300      d=>domain_glo(ind_d) 
    186301      DO j=d%jj_begin-1,d%jj_end+1 
    187302        DO i=d%ii_begin-1,d%ii_end+1 
    … …  
    203318    TYPE(t_domain),POINTER :: d  
    204319     
    205     DO ind_d=1,ndomain 
    206       d=>domain(ind_d) 
     320    DO ind_d=1,ndomain_glo 
     321      d=>domain_glo(ind_d) 
    207322      d%t_right=1 
    208323      d%t_left=-1 
    … …  
    252367  END SUBROUTINE set_neighbour_indice 
    253368       
    254        
     369  SUBROUTINE assign_domain 
     370  USE mpipara 
     371  IMPLICIT NONE 
     372    INTEGER :: nb_domain(0:mpi_size-1) 
     373    INTEGER :: rank, ind,ind_glo 
     374     
     375    DO rank=0,mpi_size-1 
     376      nb_domain(rank)=ndomain_glo/mpi_size 
     377      IF ( rank < MOD(ndomain_glo,mpi_size) ) nb_domain(rank)=nb_domain(rank)+1 
     378    ENDDO 
     379     
     380    ndomain=nb_domain(mpi_rank) 
     381    ALLOCATE(domain(ndomain)) 
     382    ALLOCATE(domloc_glo_ind(ndomain)) 
     383     
     384    rank=0 
     385    ind=0 
     386    DO ind_glo=1,ndomain_glo 
     387      ind=ind+1 
     388      domglo_rank(ind_glo)=rank 
     389      domglo_loc_ind(ind_glo)=ind 
     390      IF (rank==mpi_rank) THEN  
     391        CALL copy_domain(domain_glo(ind_glo),domain(ind)) 
     392        domloc_glo_ind(ind)=ind_glo 
     393      ENDIF 
     394       
     395      IF (ind==nb_domain(rank)) THEN 
     396        rank=rank+1 
     397        ind=0 
     398      ENDIF 
     399    ENDDO 
     400     
     401  END SUBROUTINE  assign_domain       
     402     
    255403           
    256404  SUBROUTINE compute_domain 
    257405  IMPLICIT NONE 
    258    
     406    CALL init_domain_param 
    259407    CALL create_domain 
    260408    CALL assign_cell 
    261409    CALL compute_boundary 
    262410    CALL set_neighbour_indice 
     411    CALL assign_domain 
    263412       
    264413  END SUBROUTINE compute_domain 

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

    r20 r30  
    11MODULE domain_param 
    22 
    3 INTEGER, PARAMETER :: nsplit_i=1 
    4 INTEGER, PARAMETER :: nsplit_j=1 
    5 INTEGER, PARAMETER :: halo=2 
     3INTEGER :: nsplit_i 
     4INTEGER :: nsplit_j 
     5INTEGER :: halo=2 
    66 
     7INTEGER, PARAMETER :: default_nsplit_i=1 
     8INTEGER, PARAMETER :: default_nsplit_j=1 
    79 
     10CONTAINS 
     11   
     12  SUBROUTINE init_domain_param 
     13  USE ioipsl 
     14  IMPLICIT NONE 
     15    nsplit_i=default_nsplit_i 
     16    nsplit_j=default_nsplit_j 
     17    CALL getin('nsplit_i',nsplit_i) 
     18    CALL getin('nsplit_j',nsplit_j) 
     19  END SUBROUTINE init_domain_param 
     20   
    821END MODULE domain_param 

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

    r20 r30  
    1919  REAL(rstd), PARAMETER :: rt=radius*0.5 
    2020  REAL(rstd), PARAMETER :: zc=5000.0 
    21   REAL(rstd), PARAMETER :: ps0=100000.0 
    22   REAL(rstd), PARAMETER :: T0=300 
    23   REAL(rstd), PARAMETER :: R_d=287.0    
    2421 
    2522  PUBLIC etat0 
    … …  
    7875  REAL(rstd) :: X2(3),X1(3) 
    7976  INTEGER :: i,j,n,l 
    80   INTEGER :: testcase 
    81  
    82     u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ps0 
    83          
    84     DO l=1, llm+1 
    85       pr = ap(l) + bp(l)*ps0 
    86       zr(l) = -R_d*T0/g*log(pr/ps0)  
    87     ENDDO     
    88  
    89     DO l=1, llm  
    90       zrl(l) = 0.5*(zr(l) + zr(l+1))  
    91     END DO 
    92  
    93     testcase=5 
    94     CALL getin('ncar_test_case',testcase) 
    95   
    96     SELECT CASE(testcase) 
    97 !--------------------------------------------- SINGLE COSINE BELL 
    98       CASE(0) 
    99         CALL cosine_bell_1(q)  
    100          
    101       CASE(1) 
    102         CALL cosine_bell_2(q)  
    103          
    104       CASE(2) 
    105         CALL cosine_bell_2(q)  
    106         DO l=1,llm 
    107           q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l) 
    108         END DO   
    109          
    110       CASE(3) 
    111         CALL slotted_cylinders(q)        
    112  
    113       CASE(4)  
    114         CALL cosine_bell_2(qxt1) 
    115         DO l = 1,llm 
    116           q(:,l) = 0.9 - 0.8*qxt1(:,l)*qxt1(:,l)   
    117         END DO  
    118         q = q + qxt1  
    119         CALL slotted_cylinders(qxt1) 
    120         q = q + qxt1  
    121         q = 1. - q*0.3   
    122           
    123       CASE(5)  ! hadley like meridional circulation  
    124         CALL hadleyq(q)  
    125          
    126       CASE DEFAULT 
    127         PRINT*,"no such initial profile" 
    128         STOP 
    129          
    130       END SELECT 
    131   
    132   CONTAINS 
     77  CHARACTER(len=255) :: ncar_adv_shape 
     78 
     79  u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ncar_p0 
     80   
     81  DO l=1, llm+1 
     82     pr = ap(l) + bp(l)*ncar_p0 
     83     zr(l) = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0)  
     84  ENDDO 
     85   
     86  DO l=1, llm  
     87     zrl(l) = 0.5*(zr(l) + zr(l+1))  
     88  END DO 
     89   
     90  ncar_adv_shape='cos_bell' 
     91  CALL getin('ncar_adv_shape',ncar_adv_shape) 
     92   
     93  SELECT CASE(TRIM(ncar_adv_shape)) 
     94     !--------------------------------------------- SINGLE COSINE BELL 
     95  CASE('const') 
     96     q=1 
     97  CASE('cos_bell') 
     98     CALL cosine_bell_1(q)  
     99      
     100  CASE('slotted_cyl') 
     101     CALL slotted_cylinders(q)   
     102      
     103  CASE('dbl_cos_bell_q1') 
     104     CALL cosine_bell_2(q)  
     105      
     106  CASE('dbl_cos_bell_q2') 
     107     CALL cosine_bell_2(q)  
     108     DO l=1,llm 
     109        q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l) 
     110     END DO 
     111      
     112  CASE('complement') 
     113     ! tracer such that, in combination with the other tracer fields  
     114     ! with weight (3/10), the sum is equal to one 
     115     CALL cosine_bell_2(qxt1) 
     116     DO l = 1,llm 
     117        q(:,l) = 0.9 - 0.8*qxt1(:,l)*qxt1(:,l)   
     118     END DO 
     119     q = q + qxt1  
     120     CALL slotted_cylinders(qxt1) 
     121     q = q + qxt1  
     122     q = 1. - q*0.3   
     123      
     124  CASE('hadley')  ! hadley like meridional circulation  
     125     CALL hadleyq(q)  
     126      
     127  CASE DEFAULT 
     128     PRINT *, 'Bad selector for variable ncar_adv_shape : <', TRIM(ncar_adv_shape),   & 
     129          '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', & 
     130          '<dbl_cos_bell_q2>, <complement>, <hadley>' 
     131     STOP 
     132      
     133  END SELECT 
     134       
     135CONTAINS 
    133136 
    134137!====================================================================== 

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

    r20 r30  
    2626    INTEGER :: field_type 
    2727    INTEGER :: data_type  
     28    INTEGER :: dim3 
     29    INTEGER :: dim4 
    2830  END TYPE t_field    
    2931 
    … …  
    5961      IF (PRESENT(dim2)) THEN 
    6062        field(ind)%ndim=4  
     63        field(ind)%dim4=dim2  
     64        field(ind)%dim3=dim1 
    6165      ELSE IF (PRESENT(dim1)) THEN 
    6266        field(ind)%ndim=3 
     67        field(ind)%dim3=dim1 
    6368      ELSE 
    6469        field(ind)%ndim=2 
    … …  
    96101    
    97102  END SUBROUTINE allocate_field 
    98    
     103 
     104  SUBROUTINE allocate_field_glo(field,field_type,data_type,dim1,dim2) 
     105  USE domain_mod 
     106  IMPLICIT NONE 
     107    TYPE(t_field),POINTER :: field(:) 
     108    INTEGER,INTENT(IN) :: field_type 
     109    INTEGER,INTENT(IN) :: data_type 
     110    INTEGER,OPTIONAL :: dim1,dim2 
     111    INTEGER :: ind 
     112    INTEGER :: ii_size,jj_size 
     113 
     114    ALLOCATE(field(ndomain_glo))     
     115 
     116    DO ind=1,ndomain_glo 
     117   
     118      IF (PRESENT(dim2)) THEN 
     119        field(ind)%ndim=4  
     120        field(ind)%dim4=dim2  
     121        field(ind)%dim3=dim1  
     122      ELSE IF (PRESENT(dim1)) THEN 
     123        field(ind)%ndim=3 
     124        field(ind)%dim3=dim1  
     125      ELSE 
     126        field(ind)%ndim=2 
     127      ENDIF 
     128     
     129     
     130      field(ind)%data_type=data_type 
     131      field(ind)%field_type=field_type 
     132     
     133      IF (field_type==field_T) THEN  
     134        jj_size=domain_glo(ind)%jjm 
     135      ELSE IF (field_type==field_U) THEN  
     136        jj_size=3*domain_glo(ind)%jjm 
     137      ELSE IF (field_type==field_Z) THEN  
     138        jj_size=2*domain_glo(ind)%jjm 
     139      ENDIF 
     140       
     141      ii_size=domain_glo(ind)%iim 
     142         
     143      IF (field(ind)%ndim==4) THEN 
     144        IF (data_type==type_integer) ALLOCATE(field(ind)%ival4d(ii_size*jj_size,dim1,dim2)) 
     145        IF (data_type==type_real)    ALLOCATE(field(ind)%rval4d(ii_size*jj_size,dim1,dim2)) 
     146        IF (data_type==type_logical) ALLOCATE(field(ind)%lval4d(ii_size*jj_size,dim1,dim2)) 
     147      ELSE IF (field(ind)%ndim==3) THEN 
     148        IF (data_type==type_integer) ALLOCATE(field(ind)%ival3d(ii_size*jj_size,dim1)) 
     149        IF (data_type==type_real)    ALLOCATE(field(ind)%rval3d(ii_size*jj_size,dim1)) 
     150        IF (data_type==type_logical) ALLOCATE(field(ind)%lval3d(ii_size*jj_size,dim1)) 
     151      ELSE IF (field(ind)%ndim==2) THEN 
     152        IF (data_type==type_integer) ALLOCATE(field(ind)%ival2d(ii_size*jj_size)) 
     153        IF (data_type==type_real)    ALLOCATE(field(ind)%rval2d(ii_size*jj_size)) 
     154        IF (data_type==type_logical) ALLOCATE(field(ind)%lval2d(ii_size*jj_size)) 
     155      ENDIF 
     156       
     157   ENDDO 
     158   
     159  END SUBROUTINE allocate_field_glo 
     160 
     161  SUBROUTINE deallocate_field(field) 
     162  USE domain_mod 
     163  IMPLICIT NONE 
     164    TYPE(t_field),POINTER :: field(:) 
     165    INTEGER :: data_type 
     166    INTEGER :: ind 
     167 
     168    DO ind=1,ndomain 
     169 
     170      data_type=field(ind)%data_type 
     171         
     172      IF (field(ind)%ndim==4) THEN 
     173        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival4d) 
     174        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval4d) 
     175        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval4d) 
     176      ELSE IF (field(ind)%ndim==3) THEN 
     177        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival3d) 
     178        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval3d) 
     179        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval3d) 
     180      ELSE IF (field(ind)%ndim==2) THEN 
     181        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival2d) 
     182        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval2d) 
     183        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval2d) 
     184      ENDIF 
     185       
     186   ENDDO 
     187   DEALLOCATE(field) 
     188        
     189  END SUBROUTINE deallocate_field 
     190 
     191  SUBROUTINE deallocate_field_glo(field) 
     192  USE domain_mod 
     193  IMPLICIT NONE 
     194    TYPE(t_field),POINTER :: field(:) 
     195    INTEGER :: data_type 
     196    INTEGER :: ind 
     197 
     198    DO ind=1,ndomain_glo 
     199 
     200      data_type=field(ind)%data_type 
     201         
     202      IF (field(ind)%ndim==4) THEN 
     203        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival4d) 
     204        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval4d) 
     205        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval4d) 
     206      ELSE IF (field(ind)%ndim==3) THEN 
     207        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival3d) 
     208        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval3d) 
     209        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval3d) 
     210      ELSE IF (field(ind)%ndim==2) THEN 
     211        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival2d) 
     212        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval2d) 
     213        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval2d) 
     214      ENDIF 
     215       
     216   ENDDO 
     217   DEALLOCATE(field) 
     218        
     219  END SUBROUTINE deallocate_field_glo 
     220     
    99221  SUBROUTINE getval_r2d(field_pt,field) 
    100222  IMPLICIT NONE   

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

    r20 r30  
    257257    REAL(rstd) :: surf_v(6) 
    258258    REAL(rstd) :: vect(3,6) 
    259     REAL(rstd) :: centr(6) 
     259    REAL(rstd) :: centr(3) 
    260260    REAL(rstd) :: vet(3),vep(3) 
    261261    INTEGER :: ind,i,j,k,n 

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

    r20 r30  
    1919       
    2020      CASE ('ncar') 
    21         CALL init_guided_ncar(dt) 
     21        CALL init_guided_ncar 
    2222         
    2323      CASE DEFAULT 
    … …  
    3030 
    3131   
    32   SUBROUTINE guided(it, f_ps, f_theta_rhodz, f_u, f_q) 
     32  SUBROUTINE guided(tt, f_ps, f_theta_rhodz, f_u, f_q) 
    3333  USE icosa 
    3434  USE guided_ncar_mod, ONLY : guided_ncar => guided 
    3535  IMPLICIT NONE 
    36     INTEGER, INTENT(IN)   :: it 
     36    REAL(rstd), INTENT(IN):: tt 
    3737    TYPE(t_field),POINTER :: f_ps(:) 
    3838    TYPE(t_field),POINTER :: f_phis(:) 
    … …  
    4545       
    4646      CASE ('ncar') 
    47         CALL guided_ncar(it, f_ps, f_theta_rhodz, f_u, f_q) 
     47        CALL guided_ncar(tt, f_ps, f_theta_rhodz, f_u, f_q) 
    4848    END SELECT  
    4949   

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

    r20 r30  
    33  PRIVATE 
    44   
    5   REAL(rstd),SAVE :: dt 
    6   INTEGER,SAVE    :: test 
     5  INTEGER,SAVE :: case_wind 
    76   
     7  REAL(rstd), PARAMETER :: alpha=0.0 ! tilt of solid-body rotation 
     8  REAL(rstd), PARAMETER :: tau = 12*daysec ! 12 days               ! see p. 16 
     9  REAL(rstd), PARAMETER :: w0_deform = 23000*pi/tau, b=0.2, ptop=25494.4  ! see p. 16  
     10  REAL(rstd), PARAMETER :: u0_hadley=40.,w0_hadley=0.15 ,ztop= 12000.  
     11  INTEGER, PARAMETER    :: K_hadley=5 
     12 
    813  PUBLIC init_guided, guided 
    914 
    1015CONTAINS 
    1116 
    12  
    13   SUBROUTINE init_guided(dt_in) 
    14   IMPLICIT NONE 
    15     REAL(rstd),INTENT(IN) :: dt_in 
    16      
    17     dt=dt_in 
    18     test=2 
    19      
     17  SUBROUTINE init_guided 
     18    IMPLICIT NONE 
     19    CHARACTER(LEN=255) :: wind 
     20    wind='deform' 
     21    CALL getin('ncar_adv_wind',wind) 
     22    SELECT CASE(TRIM(wind)) 
     23    CASE('solid') 
     24       case_wind=0 
     25    CASE('deform') 
     26       case_wind=1 
     27    CASE('hadley') 
     28       case_wind=2 
     29    CASE DEFAULT 
     30       PRINT*,'Bad selector for variable ncar_adv_wind : <', TRIM(wind),'> options are <solid>, <deform>, <hadley>' 
     31       END SELECT 
    2032  END SUBROUTINE init_guided 
    21  
    2233   
    23   SUBROUTINE guided(it, f_ps, f_theta_rhodz, f_u, f_q) 
     34  SUBROUTINE guided(tt, f_ps, f_theta_rhodz, f_u, f_q) 
    2435  USE icosa 
    2536    IMPLICIT NONE 
    26     INTEGER, INTENT(IN)   :: it 
     37    REAL(rstd), INTENT(IN):: tt 
    2738    TYPE(t_field),POINTER :: f_ps(:) 
    2839    TYPE(t_field),POINTER :: f_phis(:) 
    … …  
    3849      CALL swap_geometry(ind)  
    3950      ue = f_u(ind)   
    40       CALL wind_profile(it,ue) 
     51      CALL wind_profile(tt,ue) 
    4152    END DO     
    4253 
    … …  
    4455   
    4556 
    46   SUBROUTINE wind_profile(it,ue) 
    47   USE icosa 
    48   IMPLICIT NONE 
    49     INTEGER,INTENT(IN) :: it  
     57  SUBROUTINE wind_profile(tt,ue) 
     58    USE icosa 
     59    USE disvert_mod 
     60    IMPLICIT NONE 
     61    REAL(rstd),INTENT(IN)  :: tt ! current time  
    5062    REAL(rstd),INTENT(OUT) :: ue(iim*3*jjm,llm) 
    5163    REAL(rstd) :: lon, lat 
    5264    REAL(rstd) :: nx(3),n_norm,Velocity(3,llm) 
    53     REAL(rstd), PARAMETER :: lon0=3*pi/2,lat0=0.0 
    5465    REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx 
    5566    REAL(rstd) :: v1(3),v2(3),ny(3) 
    5667    INTEGER :: i,j,n,l 
    57     REAL(rstd) :: zrl(llm) 
     68    REAL(rstd) :: pitbytau,kk, pr, zr, u0, u1, v0 
    5869 
    59     CALL compute_zrl(zrl) 
     70    pitbytau = pi*tt/tau 
     71    kk = 10*radius/tau 
     72    u0 = 2*pi*radius/tau ! for solid-body rotation 
    6073!--------------------------------------------------------- 
    6174    DO l = 1,llm  
    62       DO j=jj_begin-1,jj_end+1 
    63         DO i=ii_begin-1,ii_end+1 
    64           n=(j-1)*iim+i 
    65           CALL compute_velocity(xyz_e(n+u_right,:),l,velocity(:,l)) 
    66           CALL cross_product2(xyz_v(n+z_rdown,:)/radius,xyz_v(n+z_rup,:)/radius,nx)        
    67           ue(n+u_right,l)=1e-10 
    68           n_norm=sqrt(sum(nx(:)**2)) 
    69           IF (n_norm>1e-30) THEN 
    70             nx=-nx/n_norm*ne(n,right) 
    71             ue(n+u_right,l)=sum(nx(:)*velocity(:,l)) 
    72             IF (ABS(ue(n+u_right,l))<1e-100) PRINT *,"ue(n+u_right) ==0",i,j,velocity(:,1) 
    73           ENDIF 
     75       pr = presnivs(l) 
     76       zr = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0)  ! reciprocal of (1) p. 13, isothermal atmosphere 
     77       u1 = w0_deform*radius/b/ptop*cos(2*pitbytau)*(exp((ptop-pr)/b/ptop)-exp((pr-ncar_p0)/b/ptop)) 
     78       v0 = -radius*w0_hadley*pi/(5.0*ztop)*(ncar_p0/pr)*cos(pi*zr/ztop)*cos(pitbytau) ! for Hadley cell 
     79 
     80       DO j=jj_begin-1,jj_end+1 
     81          DO i=ii_begin-1,ii_end+1 
     82             n=(j-1)*iim+i 
     83             CALL compute_velocity(xyz_e(n+u_right,:),l,velocity(:,l)) 
     84             CALL cross_product2(xyz_v(n+z_rdown,:)/radius,xyz_v(n+z_rup,:)/radius,nx)        
     85             ue(n+u_right,l)=1e-10 
     86             n_norm=sqrt(sum(nx(:)**2)) 
     87             IF (n_norm>1e-30) THEN 
     88                nx=-nx/n_norm*ne(n,right) 
     89                ue(n+u_right,l)=sum(nx(:)*velocity(:,l)) 
     90                IF (ABS(ue(n+u_right,l))<1e-100) PRINT *,"ue(n+u_right)==0",i,j,velocity(:,1) 
     91             ENDIF 
    7492              
    75           CALL compute_velocity(xyz_e(n+u_lup,:),l,velocity(:,l)) 
    76           CALL cross_product2(xyz_v(n+z_up,:)/radius,xyz_v(n+z_lup,:)/radius,nx) 
    77  
    78           ue(n+u_lup,l)=1e-10 
    79           n_norm=sqrt(sum(nx(:)**2)) 
    80           IF (n_norm>1e-30) THEN 
    81             nx=-nx/n_norm*ne(n,lup) 
    82             ue(n+u_lup,l)=sum(nx(:)*velocity(:,l)) 
    83           ENDIF 
    84         
    85           CALL compute_velocity(xyz_e(n+u_ldown,:),l,velocity(:,l)) 
    86           CALL cross_product2(xyz_v(n+z_ldown,:)/radius,xyz_v(n+z_down,:)/radius,nx) 
    87  
    88           ue(n+u_ldown,l)=1e-10 
    89           n_norm=sqrt(sum(nx(:)**2)) 
    90           IF (n_norm>1e-30) THEN 
    91             nx=-nx/n_norm*ne(n,ldown) 
    92             ue(n+u_ldown,l)=sum(nx(:)*velocity(:,l)) 
    93             IF (ABS(ue(n+u_ldown,l))<1e-100) PRINT *,"ue(n+u_ldown) ==0",i,j 
    94           ENDIF         
    95         ENDDO 
    96       ENDDO 
     93             CALL compute_velocity(xyz_e(n+u_lup,:),l,velocity(:,l)) 
     94             CALL cross_product2(xyz_v(n+z_up,:)/radius,xyz_v(n+z_lup,:)/radius,nx) 
     95              
     96             ue(n+u_lup,l)=1e-10 
     97             n_norm=sqrt(sum(nx(:)**2)) 
     98             IF (n_norm>1e-30) THEN 
     99                nx=-nx/n_norm*ne(n,lup) 
     100                ue(n+u_lup,l)=sum(nx(:)*velocity(:,l)) 
     101             ENDIF 
     102              
     103             CALL compute_velocity(xyz_e(n+u_ldown,:),l,velocity(:,l)) 
     104             CALL cross_product2(xyz_v(n+z_ldown,:)/radius,xyz_v(n+z_down,:)/radius,nx) 
     105              
     106             ue(n+u_ldown,l)=1e-10 
     107             n_norm=sqrt(sum(nx(:)**2)) 
     108             IF (n_norm>1e-30) THEN 
     109                nx=-nx/n_norm*ne(n,ldown) 
     110                ue(n+u_ldown,l)=sum(nx(:)*velocity(:,l)) 
     111                IF (ABS(ue(n+u_ldown,l))<1e-100) PRINT *,"ue(n+u_ldown)==0",i,j 
     112             ENDIF 
     113          ENDDO 
     114       ENDDO 
    97115    END DO 
    98   
     116     
    99117  CONTAINS 
    100118            
    … …  
    104122      INTEGER,INTENT(IN)::l 
    105123      REAL(rstd),INTENT(OUT) :: velocity(3) 
    106       REAL(rstd) :: u0 
    107       REAL(rstd)  :: e_lat(3), e_lon(3) 
     124      REAL(rstd) :: e_lat(3), e_lon(3) 
    108125      REAL(rstd) :: lon,lat 
    109       REAL(rstd),PARAMETER::alpha=0.0 
    110126      REAL(rstd) :: u,v 
    111 !--------------------------------- test 1 
    112       REAL(rstd),parameter:: nday=12,daysec=86400 
    113       REAL(rstd),PARAMETER :: tau = nday*daysec 
    114       REAL(rstd) :: tt, pitbytau,kk 
    115 !---------------------------------- hadley  
    116       INTEGER,PARAMETER::K_hadly=5 
    117       REAL(rstd),PARAMETER::u0_hadly=40.,w0_hadly=0.25 ,ztop= 12000.  
    118       REAL(rstd)::coff_hadly 
    119          
    120          
    121       tt = it*dt  
    122       pitbytau = pi*tt/tau  
    123       kk = 10*radius/tau  
    124        
     127               
    125128      CALL xyz2lonlat(x/radius,lon,lat) 
    126129      e_lat(1) = -cos(lon)*sin(lat) 
    … …  
    133136 
    134137      u = 0.0 ; v = 0.0 
    135       u0 = 2*pi*radius/tau 
    136138     
    137       SELECT CASE(test)  
    138         CASE(0)   
    139           u=u0*(cos(lat)*cos(alpha)+sin(lat)*sin(alpha)*cos(lon)) 
    140           v=-u0*sin(lon)*sin(alpha) 
    141         CASE(1)  
    142           lon = lon - 2*pitbytau  
    143           u= kk*sin(lon)*sin(lon)*sin(2*lat)*cos(pitbytau)+ u0*cos(lat)  
    144           v= kk*sin(2*lon)*cos(lat)*cos(pitbytau)  
    145         CASE(2) 
    146           coff_hadly = -radius*w0_hadly*pi/(5.0*ztop)   
    147           u= u0_hadly*cos(lat)  
    148           v= coff_hadly*cos(lat)*sin(5.*lat)*cos(pi*zrl(l)/ztop)*cos(pitbytau) 
    149         CASE DEFAULT  
    150           PRINT*,"not valid choice of wind"  
    151        END SELECT  
    152  
     139      SELECT CASE(case_wind)  
     140      CASE(0)  ! Solid-body rotation 
     141         u=u0*(cos(lat)*cos(alpha)+sin(lat)*sin(alpha)*cos(lon)) 
     142         v=-u0*sin(lon)*sin(alpha) 
     143      CASE(1)  ! 3D Deformational flow -  
     144         lon = lon-2*pitbytau 
     145         u = kk*sin(lon)*sin(lon)*sin(2*lat)*cos(pitbytau)+ u0*cos(lat) 
     146         u = u + u1*cos(lon)*cos(lat)**2 
     147         v = kk*sin(2*lon)*cos(lat)*cos(pitbytau)  
     148      CASE(2) ! Hadley-like flow 
     149         u = u0_hadley*cos(lat) 
     150         v = v0*cos(lat)*sin(5.*lat)    ! Eq. 37 p. 19 
     151      CASE DEFAULT  
     152         PRINT*,"not valid choice of wind"  
     153      END SELECT 
     154       
    153155      Velocity(:)=(u*e_lon(:)+v*e_lat(:)+1e-50) 
    154  
     156       
    155157    END SUBROUTINE compute_velocity 
    156158     
    157     SUBROUTINE compute_zrl(zrl) 
    158     USE disvert_mod 
    159     IMPLICIT NONE 
    160       REAL(rstd),INTENT(OUT) :: zrl(llm) 
    161       INTEGER :: l 
    162       REAL(rstd) :: zr(llm+1) 
    163       REAl(rstd) :: pr 
    164       REAL(rstd), PARAMETER :: T0=300 
    165       REAL(rstd), PARAMETER :: R_d=287.0    
    166       REAL(rstd), PARAMETER :: ps0=100000.0 
    167  
    168        
    169       DO    l    = 1, llm+1 
    170         pr = ap(l) + bp(l)*ps0 
    171         zr(l) = -R_d*T0/g*log(pr/ps0)  
    172       ENDDO     
    173  
    174       DO l = 1, llm  
    175         zrl(l) = 0.5*(zr(l) + zr(l+1))  
    176       END DO 
    177      
    178     END SUBROUTINE compute_zrl 
    179    
    180159  END SUBROUTINE  wind_profile 
    181160 

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

    r20 r30  
    22  USE icosa 
    33  USE timeloop_gcm_mod 
    4   USE transfert_mod 
    54  USE disvert_mod 
    65  USE etat0_mod 
    76  USE wind_mod 
     7  USE mpipara 
    88  IMPLICIT NONE 
    99   
    1010  TYPE(t_field),POINTER :: sum_ne(:) 
     11  TYPE(t_field),POINTER :: sum_ne_glo(:) 
    1112  REAL(rstd),POINTER :: pt_sum_ne(:) 
    1213   
    … …  
    1617  REAL(rstd) :: centr(3),dist 
    1718  
     19  CALL init_mpipara 
     20   
    1821  CALL init_grid_param 
    1922  CALL compute_metric 
    2023  CALL compute_domain 
    2124  CALL init_transfert 
     25!  CALL allocate_field(sum_ne,field_T,type_real) 
     26!  CALL allocate_field_glo(sum_ne_glo,field_T,type_real) 
     27!   
     28! DO ind=1,ndomain 
     29!   CALL swap_dimensions(ind) 
     30!   pt_sum_ne=sum_ne(ind) 
     31!   DO j=jj_begin,jj_end 
     32!     DO i=ii_begin,ii_end     
     33!       n=(j-1)*iim+i 
     34!       pt_sum_ne(n)=domloc_glo_ind(ind) 
     35!     ENDDO 
     36!   ENDDO 
     37! ENDDO 
     38 
     39! CALL WriteField("domain",sum_ne) 
     40! CALL WriteField_mpi("domain",sum_ne) 
     41! CALL transfert_request(sum_ne,req_i1) 
     42! CALL WriteField_mpi("domain",sum_ne) 
     43! CALL close_files 
     44! CALL finalize_mpipara 
     45! STOP 
     46   
    2247  CALL compute_geometry 
    2348  CALL init_disvert  
    … …  
    5277  CALL WriteField("Ai",geom%Ai) 
    5378!  CALL WriteField("sum_ne",sum_ne) 
    54    
    5579  CALL timeloop 
     80 
     81  CALL close_files 
     82  CALL finalize_mpipara 
    5683   
    5784END PROGRAM ICOSA_gcm  

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

    r20 r30  
    1515  USE transfert_mod 
    1616   
     17  ! Variables defined by run.def 
     18 
     19  REAL(rstd) :: ncar_dz, ncar_p0, ncar_T0 ! read from run.def by disvert 
     20 
    1721END MODULE icosa 

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

    r20 r30  
    22  USE icosa 
    33  USE timeloop_sw_mod 
    4   USE transfert_mod 
    54  USE dissip_sw_mod 
    65  USE disvert_mod 

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

    r20 r30  
    2525   INTEGER :: e1 
    2626   INTEGER :: e2 
     27   INTEGER :: assign_domain 
     28   INTEGER :: assign_i 
     29   INTEGER :: assign_j 
     30   INTEGER :: assign_pos 
    2731  END TYPE t_edge_glo 
    2832     

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

    r20 r30  
    9797    PRINT *,"It No :",It 
    9898 
    99     CALL guided(it,f_ps,f_theta_rhodz,f_u,f_q) 
     99    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q) 
    100100    CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du) 
    101101    CALL advect_tracer(f_ps,f_u,f_q) 
    … …  
    116116 
    117117      CASE default 
    118         PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),"> options are <euler>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>"  
     118        PRINT*,'Bad selector for variable scheme : <', TRIM(scheme),             & 
     119               ' > options are <euler>, <leapfrog>, <leapfrog_matsuno>, <adam_bashforth>' 
    119120        STOP 
    120121         

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

    r20 r30  
    1818        DO i=ii_begin,ii_end 
    1919          ij=(j-1)*iim+i 
    20           ucenter(ij,:,l)=1/Ai(ij)*(  ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_i(ij+z_rup,:))/2  -centroid(ij,:))   & 
    21                                     + ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_i(ij+z_up,:))/2-centroid(ij,:))             & 
    22                                     + ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_i(ij+z_lup,:))/2-centroid(ij,:))             & 
    23                                     + ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)*((xyz_i(ij+z_lup,:)+xyz_i(ij+z_ldown,:))/2-centroid(ij,:))       & 
    24                                     + ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)*((xyz_i(ij+z_ldown,:)+xyz_i(ij+z_down,:))/2-centroid(ij,:))   & 
    25                                     + ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown)*((xyz_i(ij+z_down,:)+xyz_i(ij+z_rdown,:))/2-centroid(ij,:)) ) 
     20          ucenter(ij,:,l)=1/Ai(ij)*                                                                                                & 
     21                        ( ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_i(ij+z_rup,:))/2  -centroid(ij,:))& 
     22                         + ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_i(ij+z_up,:))/2-centroid(ij,:))          & 
     23                         + ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_i(ij+z_lup,:))/2-centroid(ij,:))          & 
     24                         + ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)*((xyz_i(ij+z_lup,:)+xyz_i(ij+z_ldown,:))/2-centroid(ij,:))    & 
     25                         + ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)*((xyz_i(ij+z_ldown,:)+xyz_i(ij+z_down,:))/2-centroid(ij,:))& 
     26                         + ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown)*((xyz_i(ij+z_down,:)+xyz_i(ij+z_rdown,:))/2-centroid(ij,:))) 
    2627        ENDDO 
    2728      ENDDO 

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

    r20 r30  
    66    INTEGER :: size 
    77    INTEGER,POINTER :: nc_id(:) 
     8    INTEGER :: displ 
    89  END TYPE ncvar 
    910 
    … …  
    3738      enddo 
    3839    end function GetFieldIndex 
    39   
    40  
    41        
    42     subroutine WriteField_gen(name,Field,dimx,dimy,dimz) 
    43     implicit none 
    44 !    include 'netcdf.inc' 
    45       character(len=*) :: name 
    46       integer :: dimx,dimy,dimz 
    47       real,dimension(dimx,dimy,dimz) :: Field 
    48       integer,dimension(dimx*dimy*dimz) :: ndex 
    49       integer :: status 
    50       integer :: index 
    51       integer :: start(4) 
    52       integer :: count(4) 
    53        
    54             
     40     
     41    SUBROUTINE Writefield(name_in,field,nind) 
     42    USE domain_mod 
     43    USE field_mod 
     44    USE transfert_mpi_mod 
     45    USE dimensions 
     46    USE mpipara 
     47    IMPLICIT NONE   
     48     CHARACTER(LEN=*),INTENT(IN) :: name_in 
     49      TYPE(t_field),POINTER :: field(:) 
     50      INTEGER,OPTIONAL,INTENT(IN) :: nind 
     51      TYPE(t_field),POINTER :: field_glo(:) 
     52       
     53      IF (field(1)%ndim==2) THEN 
     54        CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type) 
     55      ELSE IF (field(1)%ndim==3) THEN 
     56        CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3)  
     57      ELSE IF (field(1)%ndim==4) THEN 
     58        CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3,field(1)%dim4) 
     59      ENDIF 
     60       
     61      CALL gather_field(field,field_glo) 
     62       
     63      IF (mpi_rank==0) THEN 
     64        IF (PRESENT(nind)) THEN 
     65         CALL writefield_gen(name_in,field_glo,domain_glo,nind) 
     66        ELSE 
     67         CALL writefield_gen(name_in,field_glo,domain_glo,1,ndomain_glo) 
     68        ENDIF 
     69      ENDIF 
     70       
     71      CALL deallocate_field(field_glo) 
     72       
     73   END SUBROUTINE writefield 
     74    
     75!   SUBROUTINE Writefield(name_in,field,nind) 
     76!   USE netcdf 
     77!   USE domain_mod 
     78!   use field_mod 
     79!   USE dimensions 
     80!   USE geometry 
     81!   IMPLICIT NONE 
     82!     CHARACTER(LEN=*),INTENT(IN) :: name_in 
     83!     TYPE(t_field),POINTER :: field(:) 
     84!     INTEGER,OPTIONAL,INTENT(IN) :: nind 
     85!     REAL(r8),ALLOCATABLE :: field_val2d(:) 
     86!     REAL(r8),ALLOCATABLE :: field_val3d(:,:) 
     87!     REAL(r8),ALLOCATABLE :: field_val4d(:,:,:) 
     88!     TYPE(t_domain),POINTER :: d 
     89!     INTEGER :: Index 
     90!     INTEGER :: ind,i,j,k,n,ncell,q 
     91!     INTEGER :: iie,jje,iin,jjn 
     92!     INTEGER :: status 
     93!     CHARACTER(len=255) :: name 
     94!     CHARACTER(len=255) :: str_ind 
     95!     INTEGER :: ind_b,ind_e 
     96!     INTEGER :: halo_size 
     97!     LOGICAL :: single 
     98!      
     99!      
     100!     name=TRIM(ADJUSTL(name_in)) 
     101 
     102!     IF (PRESENT(nind)) THEN 
     103!       name=TRIM(name)//"_"//TRIM(int2str(nind)) 
     104!       PRINT *,"NAME",nind,int2str(nind),name 
     105!       ind_b=nind 
     106!       ind_e=nind 
     107!       halo_size=1 
     108!       single=.TRUE. 
     109!     ELSE 
     110!       ind_b=1 
     111!       ind_e=ndomain 
     112!       halo_size=0 
     113!       single=.FALSE. 
     114!     ENDIF       
     115 
     116!     Index=GetFieldIndex(name) 
     117!     if (Index==-1) then 
     118!       call create_header(name,field,nind) 
     119!       Index=GetFieldIndex(name) 
     120!     else 
     121!       FieldIndex(Index)=FieldIndex(Index)+1. 
     122!     endif 
     123!      
     124!     IF (Field(ind_b)%field_type==field_T) THEN 
     125!       ncell=1 
     126!       DO ind=ind_b,ind_e 
     127!         d=>domain(ind) 
     128!         IF (Field(ind)%field_type/=field_T) THEN 
     129!           PRINT *,"Writefield, grille non geree" 
     130!           RETURN 
     131!         ENDIF 
     132 
     133!       n=0 
     134!       DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     135!         DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     136!           IF (d%own(i,j) .OR. single) n=n+1 
     137!         ENDDO 
     138!       ENDDO 
     139 
     140!       IF (field(ind)%ndim==2) THEN 
     141!         ALLOCATE(Field_val2d(n)) 
     142!         n=0 
     143!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     144!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     145!             k=d%iim*(j-1)+i 
     146!             IF (d%own(i,j) .OR. single) THEN 
     147!               n=n+1 
     148!               Field_val2d(n)=field(ind)%rval2d(k) 
     149!             ENDIF 
     150!           ENDDO 
     151!         ENDDO 
     152!         status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  & 
     153!                             start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
     154!         DEALLOCATE(field_val2d) 
     155!       ELSE IF (field(ind)%ndim==3) THEN 
     156!         ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2))) 
     157!         n=0 
     158!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     159!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     160!             k=d%iim*(j-1)+i 
     161!             IF (d%own(i,j) .OR. single) THEN 
     162!               n=n+1 
     163!               Field_val3d(n,:)=field(ind)%rval3d(k,:) 
     164!             ENDIF 
     165!           ENDDO 
     166!         ENDDO 
     167!          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
     168!                              count=(/n,size(field(1)%rval3d,2),1 /)) 
     169!         DEALLOCATE(field_val3d) 
     170!       ELSE IF (field(1)%ndim==4) THEN 
     171 
     172!         DO q=1,FieldVarId(index)%size 
     173!          
     174!           ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2))) 
     175!           n=0 
     176!           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     177!             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     178!               k=d%iim*(j-1)+i 
     179!               IF (d%own(i,j) .OR. single) THEN 
     180!                 n=n+1 
     181!                 Field_val3d(n,:)=field(ind)%rval4d(k,:,q) 
     182!               ENDIF 
     183!             ENDDO 
     184!           ENDDO 
     185!           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
     186!                              count=(/n,size(field(1)%rval4d,2),1 /)) 
     187!           DEALLOCATE(field_val3d) 
     188!         ENDDO          
     189!       ENDIF 
     190!        
     191!        PRINT *,NF90_STRERROR(status) 
     192!       ncell=ncell+n 
     193 
     194!    ENDDO 
     195!     
     196!    ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
     197!       ncell=1 
     198!       n=0 
     199!       DO ind=ind_b,ind_e 
     200!         d=>domain(ind) 
     201!         CALL swap_geometry(ind) 
     202!         CALL swap_dimensions(ind) 
     203!  
     204!         n=0 
     205!         DO j=jj_begin+1,jj_end 
     206!           DO i=ii_begin,ii_end-1 
     207!             n=n+1 
     208!           ENDDO 
     209!         ENDDO 
     210 
     211!         DO j=jj_begin,jj_end-1 
     212!           DO i=ii_begin+1,ii_end 
     213!             n=n+1 
     214!           ENDDO 
     215!         ENDDO 
     216 
     217!       IF (field(ind)%ndim==2) THEN 
     218!         ALLOCATE(Field_val2d(n)) 
     219 
     220!         n=0 
     221!         DO j=jj_begin+1,jj_end 
     222!           DO i=ii_begin,ii_end-1 
     223!             n=n+1 
     224!             k=iim*(j-1)+i 
     225!             Field_val2d(n)=field(ind)%rval2d(k+z_down) 
     226!           ENDDO 
     227!         ENDDO 
     228 
     229!         DO j=jj_begin,jj_end-1 
     230!           DO i=ii_begin+1,ii_end 
     231!             n=n+1 
     232!             k=iim*(j-1)+i 
     233!             Field_val2d(n)=field(ind)%rval2d(k+z_up) 
     234!           ENDDO 
     235!         ENDDO 
     236 
     237!         status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       & 
     238!                             Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
     239!         DEALLOCATE(field_val2d) 
     240 
     241!       ELSE IF (field(ind)%ndim==3) THEN 
     242!         ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2))) 
     243!         n=0 
     244!         DO j=jj_begin+1,jj_end 
     245!           DO i=ii_begin,ii_end-1 
     246!             n=n+1 
     247!             k=iim*(j-1)+i 
     248!             Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:) 
     249!           ENDDO 
     250!         ENDDO 
     251 
     252!         DO j=jj_begin,jj_end-1 
     253!           DO i=ii_begin+1,ii_end 
     254!             n=n+1 
     255!             k=iim*(j-1)+i 
     256!             Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:) 
     257!           ENDDO 
     258!         ENDDO 
     259!          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
     260!                              count=(/n,size(field(1)%rval3d,2),1 /)) 
     261!         DEALLOCATE(field_val3d) 
     262!       ELSE IF (field(1)%ndim==4) THEN 
     263 
     264!         DO q=1,FieldVarId(index)%size 
     265!           ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2))) 
     266!           n=0 
     267!           DO j=jj_begin+1,jj_end 
     268!             DO i=ii_begin,ii_end-1 
     269!               n=n+1 
     270!               k=iim*(j-1)+i 
     271!               Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q) 
     272!             ENDDO 
     273!           ENDDO 
     274 
     275!           DO j=jj_begin,jj_end-1 
     276!             DO i=ii_begin+1,ii_end 
     277!               n=n+1 
     278!               k=iim*(j-1)+i 
     279!               Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q) 
     280!             ENDDO 
     281!           ENDDO 
     282 
     283!           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,1,FieldIndex(Index) /),   & 
     284!                               count=(/n,size(field(1)%rval4d,2),1 /)) 
     285!           DEALLOCATE(field_val3d) 
     286!         ENDDO 
     287!       ENDIF 
     288!        
     289!        PRINT *,NF90_STRERROR(status) 
     290!       ncell=ncell+n 
     291 
     292!    ENDDO 
     293!     
     294!    ENDIF 
     295!    status=NF90_SYNC(FieldId(Index)) 
     296!      
     297!   END SUBROUTINE Writefield 
     298 
     299 
     300    SUBROUTINE Writefield_gen(name_in, field, domain_type, ind_b_in, ind_e_in ) 
     301    USE netcdf_mod 
     302    USE domain_mod 
     303    USE field_mod 
     304!    USE dimensions 
     305!    USE geometry 
     306    IMPLICIT NONE 
     307      CHARACTER(LEN=*),INTENT(IN) :: name_in 
     308      TYPE(t_field), POINTER :: field(:) 
     309      TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:) 
     310      INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in 
     311      INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in 
     312      REAL(r8),ALLOCATABLE :: field_val2d(:) 
     313      REAL(r8),ALLOCATABLE :: field_val3d(:,:) 
     314      REAL(r8),ALLOCATABLE :: field_val4d(:,:,:) 
     315      TYPE(t_domain),POINTER :: d 
     316      INTEGER :: Index 
     317      INTEGER :: ind,i,j,k,n,ncell,q 
     318      INTEGER :: iie,jje,iin,jjn 
     319      INTEGER :: status 
     320      CHARACTER(len=255) :: name 
     321      CHARACTER(len=255) :: str_ind 
     322      INTEGER :: ind_b,ind_e 
     323      INTEGER :: halo_size 
     324      LOGICAL :: single 
     325       
     326       
     327      name=TRIM(ADJUSTL(name_in)) 
     328 
     329      IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN 
     330        name=TRIM(name)//"_"//TRIM(int2str(ind_b)) 
     331        ind_b=ind_b_in 
     332        ind_e=ind_b_in 
     333        halo_size=1 
     334        single=.TRUE. 
     335      ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
     336        name=TRIM(name)//"_"//TRIM(int2str(ind_e)) 
     337        ind_b=ind_e_in 
     338        ind_e=ind_e_in 
     339        halo_size=1 
     340        single=.TRUE. 
     341      ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
     342        ind_b=ind_b_in 
     343        ind_e=ind_e_in 
     344        halo_size=0 
     345        single=.FALSE. 
     346      ELSE 
     347        halo_size=0 
     348        ind_b=1 
     349        ind_e=ndomain 
     350        single=.FALSE. 
     351      ENDIF       
     352       
    55353      Index=GetFieldIndex(name) 
    56354      if (Index==-1) then 
    57 !        call CreateNewField(name,dimx,dimy,dimz) 
    58         Index=GetFieldIndex(name) 
     355        call create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in) 
     356        Index=GetFieldIndex(name) 
    59357      else 
    60358        FieldIndex(Index)=FieldIndex(Index)+1. 
    61359      endif 
    62360       
    63       start(1)=1 
    64       start(2)=1 
    65       start(3)=1 
    66       start(4)=FieldIndex(Index) 
    67  
    68       count(1)=dimx 
    69       count(2)=dimy 
    70       count(3)=dimz 
    71       count(4)=1 
    72  
    73 !      status = NF_PUT_VARA_DOUBLE(FieldId(Index),FieldVarId(Index),start,count,Field) 
    74 !      status = NF_SYNC(FieldId(Index)) 
    75        
    76     end subroutine WriteField_gen 
    77  
    78  
    79     SUBROUTINE Writefield(name_in,field,nind) 
    80     USE netcdf 
     361      IF (Field(ind_b)%field_type==field_T) THEN 
     362 
     363        ncell=0 
     364        DO ind=ind_b,ind_e 
     365          d=>domain_type(ind) 
     366          DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     367            DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     368              IF (d%assign_domain(i,j)==ind .OR. single) ncell=ncell+1 
     369            ENDDO 
     370          ENDDO 
     371        ENDDO 
     372         
     373        IF (field(1)%ndim==2) THEN 
     374          ALLOCATE(Field_val2d(ncell)) 
     375          n=0 
     376          DO ind=ind_b,ind_e 
     377            d=>domain_type(ind) 
     378            DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     379              DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     380                k=d%iim*(j-1)+i 
     381                IF (d%assign_domain(i,j)==ind .OR. single) THEN 
     382                  n=n+1 
     383                  Field_val2d(n)=field(ind)%rval2d(k) 
     384                ENDIF 
     385              ENDDO 
     386            ENDDO 
     387          ENDDO 
     388          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  & 
     389                              start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /)) 
     390          DEALLOCATE(field_val2d) 
     391        ELSE IF (field(1)%ndim==3) THEN 
     392          ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2))) 
     393          n=0 
     394          DO ind=ind_b,ind_e 
     395            d=>domain_type(ind) 
     396            DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     397              DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     398                k=d%iim*(j-1)+i 
     399                IF (d%assign_domain(i,j)==ind .OR. single) THEN 
     400                  n=n+1 
     401                  Field_val3d(n,:)=field(ind)%rval3d(k,:) 
     402                ENDIF 
     403              ENDDO 
     404            ENDDO 
     405          ENDDO 
     406            
     407          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   & 
     408                               count=(/ncell,size(field(1)%rval3d,2),1 /)) 
     409          DEALLOCATE(field_val3d) 
     410        ELSE IF (field(1)%ndim==4) THEN 
     411 
     412          DO q=1,FieldVarId(index)%size 
     413           
     414            ALLOCATE(Field_val3d(ncell,size(field(1)%rval4d,2))) 
     415            n=0 
     416            DO ind=ind_b,ind_e 
     417              d=>domain_type(ind) 
     418              DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     419                DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     420                  k=d%iim*(j-1)+i 
     421                  IF (d%assign_domain(i,j)==ind .OR. single) THEN 
     422                    n=n+1 
     423                    Field_val3d(n,:)=field(ind)%rval4d(k,:,q) 
     424                  ENDIF 
     425                ENDDO 
     426              ENDDO 
     427            ENDDO 
     428             
     429            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   & 
     430                               count=(/ncell,size(field(1)%rval4d,2),1 /)) 
     431            DEALLOCATE(field_val3d) 
     432          ENDDO          
     433        ENDIF 
     434         
     435!        PRINT *,NF90_STRERROR(status) 
     436!        ncell=ncell+n 
     437 
     438!     ENDDO 
     439      
     440     ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
     441 
     442        ncell=0 
     443        DO ind=ind_b,ind_e 
     444          d=>domain_type(ind) 
     445   
     446          DO j=d%jj_begin+1,d%jj_end 
     447            DO i=d%ii_begin,d%ii_end-1 
     448              ncell=ncell+1 
     449            ENDDO 
     450          ENDDO 
     451 
     452          DO j=d%jj_begin,d%jj_end-1 
     453            DO i=d%ii_begin+1,d%ii_end 
     454              ncell=ncell+1 
     455            ENDDO 
     456          ENDDO 
     457        ENDDO 
     458 
     459        IF (field(1)%ndim==2) THEN 
     460          ALLOCATE(Field_val2d(ncell)) 
     461 
     462          n=0 
     463           
     464          DO ind=ind_b,ind_e 
     465            d=>domain_type(ind) 
     466            DO j=d%jj_begin+1,d%jj_end 
     467              DO i=d%ii_begin,d%ii_end-1 
     468                n=n+1 
     469                k=d%iim*(j-1)+i 
     470                Field_val2d(n)=field(ind)%rval2d(k+d%z_down) 
     471              ENDDO 
     472            ENDDO 
     473 
     474            DO j=d%jj_begin,d%jj_end-1 
     475              DO i=d%ii_begin+1,d%ii_end 
     476                n=n+1 
     477                k=d%iim*(j-1)+i 
     478                Field_val2d(n)=field(ind)%rval2d(k+d%z_up) 
     479              ENDDO 
     480            ENDDO 
     481          ENDDO 
     482           
     483          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       & 
     484                              Field_val2d,start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /)) 
     485          DEALLOCATE(field_val2d) 
     486 
     487        ELSE IF (field(1)%ndim==3) THEN 
     488          ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2))) 
     489          n=0 
     490          DO ind=ind_b,ind_e 
     491            d=>domain_type(ind) 
     492            DO j=d%jj_begin+1,d%jj_end 
     493              DO i=d%ii_begin,d%ii_end-1 
     494                n=n+1 
     495                k=d%iim*(j-1)+i 
     496                Field_val3d(n,:)=field(ind)%rval3d(k+d%z_down,:) 
     497              ENDDO 
     498            ENDDO 
     499 
     500            DO j=d%jj_begin,d%jj_end-1 
     501              DO i=d%ii_begin+1,d%ii_end 
     502                n=n+1 
     503                k=d%iim*(j-1)+i 
     504                Field_val3d(n,:)=field(ind)%rval3d(k+d%z_up,:) 
     505              ENDDO 
     506            ENDDO 
     507          ENDDO 
     508           
     509          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   & 
     510                              count=(/ncell,size(field(1)%rval3d,2),1 /)) 
     511          DEALLOCATE(field_val3d) 
     512         
     513        ELSE IF (field(1)%ndim==4) THEN 
     514 
     515          DO q=1,FieldVarId(index)%size 
     516            ALLOCATE(Field_val3d(ncell,size(field(1)%rval4d,2))) 
     517            n=0 
     518            DO ind=ind_b,ind_e 
     519              d=>domain_type(ind) 
     520              DO j=d%jj_begin+1,d%jj_end 
     521                DO i=d%ii_begin,d%ii_end-1 
     522                  n=n+1 
     523                  k=d%iim*(j-1)+i 
     524                  Field_val3d(n,:)=field(ind)%rval4d(k+d%z_down,:,q) 
     525                ENDDO 
     526              ENDDO 
     527 
     528              DO j=d%jj_begin,d%jj_end-1 
     529                DO i=d%ii_begin+1,d%ii_end 
     530                  n=n+1 
     531                  k=d%iim*(j-1)+i 
     532                  Field_val3d(n,:)=field(ind)%rval4d(k+d%z_up,:,q) 
     533                ENDDO 
     534              ENDDO 
     535            ENDDO 
     536 
     537            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,1,FieldIndex(Index) /),   & 
     538                                count=(/ncell,size(field(1)%rval4d,2),1 /)) 
     539            DEALLOCATE(field_val3d) 
     540          ENDDO 
     541        ENDIF 
     542         
     543!        PRINT *,NF90_STRERROR(status) 
     544!        ncell=ncell+n 
     545! 
     546!     ENDDO 
     547      
     548     ENDIF 
     549     status=NF90_SYNC(FieldId(Index)) 
     550       
     551    END SUBROUTINE Writefield_gen 
     552 
     553       
     554 
     555    SUBROUTINE Writefield_mpi(name_in,field,nind) 
     556    USE netcdf_mod 
    81557    USE domain_mod 
    82558    use field_mod 
    … …  
    92568      TYPE(t_domain),POINTER :: d 
    93569      INTEGER :: Index 
    94       INTEGER :: ind,i,j,k,n,ncell,q 
     570      INTEGER :: ind,i,j,l,k,n,ncell,q 
    95571      INTEGER :: iie,jje,iin,jjn 
    96572      INTEGER :: status 
    … …  
    100576      INTEGER :: halo_size 
    101577      LOGICAL :: single 
     578      INTEGER :: displ 
    102579       
    103580       
    … …  
    120597      Index=GetFieldIndex(name) 
    121598      if (Index==-1) then 
    122         call create_header(name,field,nind) 
     599        call create_header_mpi(name,field,nind) 
    123600        Index=GetFieldIndex(name) 
    124601      else 
    … …  
    142619        ENDDO 
    143620 
     621        displ=FieldVarId(index)%displ 
     622         
    144623        IF (field(ind)%ndim==2) THEN 
    145624          ALLOCATE(Field_val2d(n)) 
    … …  
    154633            ENDDO 
    155634          ENDDO 
    156           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
     635          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  & 
     636                              start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /)) 
    157637          DEALLOCATE(field_val2d) 
    158638        ELSE IF (field(ind)%ndim==3) THEN 
    … …  
    168648            ENDDO 
    169649          ENDDO 
    170            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
    171                                count=(/n,size(field(1)%rval3d,2),1 /)) 
     650!           DO l=1,size(field(ind)%rval3d,2) 
     651            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ displ+ncell,1,FieldIndex(Index) /),   & 
     652                               count=(/n,size(field(ind)%rval3d,2),1 /)) 
     653!           ENDDO 
    172654          DEALLOCATE(field_val3d) 
    173655        ELSE IF (field(1)%ndim==4) THEN 
    … …  
    185667                ENDIF 
    186668              ENDDO 
    187             ENDDO 
    188             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
    189                                count=(/n,size(field(1)%rval4d,2),1 /)) 
     669             ENDDO 
     670!            DO l=1,size(field(ind)%rval4d,2) 
     671            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d(:,l),start=(/ displ+ncell,l,FieldIndex(Index) /),   & 
     672                               count=(/n,size(field(ind)%rval4d,2),1 /)) 
     673!            ENDDO 
    190674            DEALLOCATE(field_val3d) 
    191675          ENDDO          
    … …  
    218702          ENDDO 
    219703 
     704        displ=FieldVarId(index)%displ 
     705 
    220706        IF (field(ind)%ndim==2) THEN 
    221707          ALLOCATE(Field_val2d(n)) 
    … …  
    238724          ENDDO 
    239725 
    240           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /)) 
     726          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       & 
     727                              Field_val2d,start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /)) 
    241728          DEALLOCATE(field_val2d) 
    242729 
    … …  
    259746            ENDDO 
    260747          ENDDO 
    261            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   & 
    262                                count=(/n,size(field(1)%rval3d,2),1 /)) 
     748!           DO l=1,size(field(ind)%rval3d,2) 
     749           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ displ+ncell,1,FieldIndex(Index) /),   & 
     750                               count=(/n,size(field(ind)%rval3d,2),1 /)) 
     751!           ENDDO 
    263752          DEALLOCATE(field_val3d) 
    264753        ELSE IF (field(1)%ndim==4) THEN 
    … …  
    283772            ENDDO 
    284773 
    285             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,1,FieldIndex(Index) /),   & 
    286                                 count=(/n,size(field(1)%rval4d,2),1 /)) 
     774!            DO l=1,size(field(ind)%rval4d,2) 
     775 
     776            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ displ+ncell,1,FieldIndex(Index) /),   & 
     777                                count=(/n,size(field(ind)%rval4d,2),1 /)) 
     778!            ENDDO 
    287779            DEALLOCATE(field_val3d) 
    288780          ENDDO 
    … …  
    297789     status=NF90_SYNC(FieldId(Index)) 
    298790       
    299     END SUBROUTINE Writefield 
    300        
     791    END SUBROUTINE Writefield_mpi 
     792     
     793     
    301794            
    302     SUBROUTINE Create_header(name,field,nind) 
    303     USE netcdf 
     795!   SUBROUTINE Create_header(name,field,nind) 
     796!   USE netcdf 
     797!   USE field_mod 
     798!   USE domain_mod 
     799!   USE spherical_geom_mod 
     800!   USE dimensions 
     801!   USE geometry 
     802!   IMPLICIT NONE 
     803!     CHARACTER(LEN=*) :: name 
     804!     TYPE(t_field),POINTER :: field(:) 
     805!     INTEGER,OPTIONAL,INTENT(IN) :: nind 
     806!     INTEGER :: ncell 
     807!     INTEGER :: nvert 
     808!     REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:) 
     809!     TYPE(t_domain),POINTER :: d 
     810!     INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId 
     811!     INTEGER :: dim3id,dim4id 
     812!     INTEGER :: status 
     813!     INTEGER :: ind,i,j,k,n,q 
     814!     INTEGER :: iie,jje,iin,jjn 
     815!     INTEGER :: ind_b,ind_e 
     816!     INTEGER :: halo_size 
     817!     LOGICAL :: single  
     818!     INTEGER :: nij 
     819!          
     820!     NbField=NbField+1 
     821!     FieldName(NbField)=TRIM(ADJUSTL(name)) 
     822!     FieldIndex(NbField)=1 
     823!      
     824!     IF (PRESENT(nind)) THEN 
     825!       ind_b=nind 
     826!       ind_e=nind 
     827!       halo_size=1 
     828!       single=.TRUE. 
     829!     ELSE 
     830!       ind_b=1 
     831!       ind_e=ndomain 
     832!       halo_size=0 
     833!       single=.FALSE. 
     834!     ENDIF 
     835!      
     836!     ncell=0 
     837!      
     838!     IF (Field(ind_b)%field_type==field_T) THEN 
     839!       nvert=6 
     840!        
     841!       DO ind=ind_b,ind_e 
     842!         d=>domain(ind) 
     843!        
     844!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     845!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     846!             IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1 
     847!           ENDDO 
     848!         ENDDO 
     849 
     850!       END DO 
     851!      
     852!       status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     853!       FieldId(NbField)=ncid 
     854!       status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     855!       status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
     856 
     857!       IF (Field(ind_b)%ndim==2)  THEN 
     858!         FieldVarId(NbField)%size=1 
     859!         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     860!       ELSE IF (Field(ind_b)%ndim==3)  THEN 
     861!         FieldVarId(NbField)%size=1 
     862!         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     863!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
     864!       ELSE IF (Field(1)%ndim==4) THEN 
     865!         FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
     866!         ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
     867!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
     868!          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) 
     869!       ENDIF 
     870!      
     871!       status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
     872!      
     873!       status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
     874!       status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
     875!       status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
     876!       status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
     877!       status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
     878!       status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
     879!       status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
     880!       status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
     881!       status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
     882!       status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
     883 
     884!       IF (Field(ind_b)%ndim==2) THEN 
     885!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
     886!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     887!       ELSE IF (Field(ind_b)%ndim==3) THEN 
     888!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
     889!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     890!       ELSE IF (Field(ind_b)%ndim==4) THEN 
     891!         DO i=1,FieldVarId(NbField)%size 
     892!           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),  & 
     893!                                 FieldVarId(NbField)%nc_id(i)) 
     894!           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat") 
     895!         ENDDO         
     896!       ENDIF  
     897!  
     898!      
     899!       status = NF90_ENDDEF(ncid)       
     900 
     901!       ncell=1 
     902!       DO ind=ind_b,ind_e 
     903!         d=>domain(ind) 
     904!  
     905!         n=0 
     906!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     907!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     908!             IF (single .OR. d%own(i,j)) n=n+1 
     909!           ENDDO 
     910!         ENDDO 
     911!          
     912!        ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) 
     913!          
     914!         n=0   
     915!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     916!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     917!               IF (d%own(i,j) .OR. single) THEN 
     918!               n=n+1 
     919!               CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n)) 
     920!               lon(n)=lon(n)*180/Pi 
     921!               lat(n)=lat(n)*180/Pi 
     922!               DO k=0,5 
     923!                 CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n)) 
     924!                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     925!                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     926!               ENDDO 
     927!             ENDIF 
     928!           ENDDO 
     929!         ENDDO 
     930!         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) 
     931!         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) 
     932!         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     933!         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     934!  
     935!         ncell=ncell+n 
     936!         DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
     937!     END DO  
     938 
     939!   ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
     940!       nvert=3 
     941!       DO ind=ind_b,ind_e 
     942!         d=>domain(ind) 
     943!        
     944!         DO j=d%jj_begin+1,d%jj_end 
     945!           DO i=d%ii_begin,d%ii_end-1 
     946!             ncell=ncell+1 
     947!           ENDDO 
     948!         ENDDO 
     949 
     950!         DO j=d%jj_begin,d%jj_end-1 
     951!           DO i=d%ii_begin+1,d%ii_end 
     952!             ncell=ncell+1 
     953!           ENDDO 
     954!         ENDDO 
     955 
     956!       END DO 
     957!      
     958!       status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     959!       FieldId(NbField)=ncid 
     960!       status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     961!       status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
     962 
     963!       IF (Field(ind_b)%ndim==2)  THEN 
     964!         FieldVarId(NbField)%size=1 
     965!         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     966!       ELSE IF (Field(ind_b)%ndim==3)  THEN 
     967!         FieldVarId(NbField)%size=1 
     968!         ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     969!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
     970!       ELSE IF (Field(1)%ndim==4) THEN 
     971!         FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
     972!         ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
     973!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
     974!       ENDIF 
     975 
     976 
     977!      
     978!       status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
     979!      
     980!       status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
     981!       status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
     982!       status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
     983!       status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
     984!       status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
     985!       status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
     986!       status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
     987!       status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
     988!       status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
     989!       status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
     990 
     991 
     992!       IF (Field(ind_b)%ndim==2) THEN 
     993!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
     994!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     995!       ELSE IF (Field(ind_b)%ndim==3) THEN 
     996!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
     997!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     998!       ELSE IF (Field(ind_b)%ndim==4) THEN 
     999!         DO q=1,FieldVarId(NbField)%size 
     1000!           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE,             & 
     1001!                                 (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) 
     1002!           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat") 
     1003!         ENDDO         
     1004!       ENDIF  
     1005!        
     1006!       status = NF90_ENDDEF(ncid)       
     1007 
     1008!       ncell=1 
     1009!       DO ind=ind_b,ind_e 
     1010!         d=>domain(ind) 
     1011!         CALL swap_geometry(ind) 
     1012!         CALL swap_dimensions(ind) 
     1013!  
     1014!         n=0 
     1015!         DO j=jj_begin+1,jj_end 
     1016!           DO i=ii_begin,ii_end-1 
     1017!             n=n+1 
     1018!           ENDDO 
     1019!         ENDDO 
     1020 
     1021!         DO j=jj_begin,jj_end-1 
     1022!           DO i=ii_begin+1,ii_end 
     1023!             n=n+1 
     1024!           ENDDO 
     1025!         ENDDO 
     1026 
     1027!        ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) 
     1028!          
     1029!         n=0   
     1030!       
     1031!         DO j=jj_begin+1,jj_end 
     1032!           DO i=ii_begin,ii_end-1 
     1033!             nij=(j-1)*iim+i 
     1034!             n=n+1 
     1035!             CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n)) 
     1036!             lon(n)=lon(n)*180/Pi 
     1037!!             IF (lon(n)<0) lon(n)=lon(n)+360 
     1038!             lat(n)=lat(n)*180/Pi 
     1039!             CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
     1040!             CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
     1041!             CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
     1042 
     1043!             DO k=0,2 
     1044!               bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     1045!               bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     1046!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
     1047!             ENDDO 
     1048!           ENDDO 
     1049!         ENDDO 
     1050 
     1051!         DO j=jj_begin,jj_end-1 
     1052!           DO i=ii_begin+1,ii_end 
     1053!             nij=(j-1)*iim+i 
     1054!             n=n+1 
     1055!             CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n)) 
     1056!             lon(n)=lon(n)*180/Pi 
     1057!              IF (lon(n)<0) lon(n)=lon(n)+360 
     1058!             lat(n)=lat(n)*180/Pi 
     1059!             CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
     1060!             CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
     1061!             CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
     1062 
     1063!             DO k=0,2 
     1064!               bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     1065!               bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     1066!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
     1067!             ENDDO 
     1068!           ENDDO 
     1069!         ENDDO 
     1070!          
     1071!          
     1072!         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) 
     1073!         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) 
     1074!         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     1075!         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     1076!  
     1077!         ncell=ncell+n 
     1078!         DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
     1079!     END DO           
     1080!   ENDIF 
     1081 
     1082 
     1083!      
     1084!      status = NF90_CLOSE(ncid) 
     1085 
     1086!   END SUBROUTINE Create_Header 
     1087 
     1088 
     1089 
     1090    
     1091    SUBROUTINE Create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in) 
     1092    USE netcdf_mod 
     1093    USE field_mod 
     1094    USE domain_mod 
     1095    USE metric 
     1096    USE spherical_geom_mod 
     1097!    USE dimensions 
     1098!    USE geometry 
     1099    IMPLICIT NONE 
     1100      CHARACTER(LEN=*),INTENT(IN) :: name_in 
     1101      TYPE(t_field),POINTER :: field(:) 
     1102      TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:) 
     1103      INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in 
     1104      INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in 
     1105      INTEGER :: ncell 
     1106      INTEGER :: nvert 
     1107      REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:) 
     1108      TYPE(t_domain),POINTER :: d 
     1109      INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId 
     1110      INTEGER :: dim3id,dim4id 
     1111      INTEGER :: status 
     1112      INTEGER :: ind,i,j,k,n,q 
     1113      INTEGER :: iie,jje,iin,jjn 
     1114      INTEGER :: ind_b,ind_e 
     1115      INTEGER :: halo_size 
     1116      LOGICAL :: single  
     1117      INTEGER :: nij 
     1118      CHARACTER(LEN=255) :: name 
     1119       
     1120             
     1121      name=TRIM(ADJUSTL(name_in)) 
     1122 
     1123      IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN 
     1124        name=TRIM(name)//"_"//TRIM(int2str(ind_b)) 
     1125        ind_b=ind_b_in 
     1126        ind_e=ind_b_in 
     1127        halo_size=1 
     1128        single=.TRUE. 
     1129      ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
     1130        name=TRIM(name)//"_"//TRIM(int2str(ind_e)) 
     1131        ind_b=ind_e_in 
     1132        ind_e=ind_e_in 
     1133        halo_size=1 
     1134        single=.TRUE. 
     1135      ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN 
     1136        ind_b=ind_b_in 
     1137        ind_e=ind_e_in 
     1138        halo_size=0 
     1139        single=.FALSE. 
     1140      ELSE 
     1141        halo_size=0 
     1142        ind_b=1 
     1143        ind_e=ndomain 
     1144        single=.FALSE. 
     1145      ENDIF       
     1146                 
     1147      NbField=NbField+1 
     1148      FieldName(NbField)=TRIM(ADJUSTL(name)) 
     1149      FieldIndex(NbField)=1 
     1150       
     1151      ncell=0 
     1152       
     1153      IF (Field(ind_b)%field_type==field_T) THEN 
     1154        nvert=6 
     1155         
     1156        DO ind=ind_b,ind_e 
     1157          d=>domain_type(ind) 
     1158         
     1159          DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     1160            DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     1161              IF (single .OR. d%assign_domain(i,j)==ind) ncell=ncell+1 
     1162            ENDDO 
     1163          ENDDO 
     1164 
     1165        END DO 
     1166       
     1167        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     1168        FieldId(NbField)=ncid 
     1169        status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     1170        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
     1171 
     1172        IF (Field(ind_b)%ndim==2)  THEN 
     1173          FieldVarId(NbField)%size=1 
     1174          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     1175        ELSE IF (Field(ind_b)%ndim==3)  THEN 
     1176          FieldVarId(NbField)%size=1 
     1177          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     1178          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
     1179        ELSE IF (Field(1)%ndim==4) THEN 
     1180          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
     1181          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
     1182          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
     1183!          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id) 
     1184        ENDIF 
     1185       
     1186        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
     1187       
     1188        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
     1189        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
     1190        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
     1191        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
     1192        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
     1193        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
     1194        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
     1195        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
     1196        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
     1197        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
     1198 
     1199        IF (Field(ind_b)%ndim==2) THEN 
     1200          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
     1201          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1202        ELSE IF (Field(ind_b)%ndim==3) THEN 
     1203          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
     1204          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1205        ELSE IF (Field(ind_b)%ndim==4) THEN 
     1206          DO i=1,FieldVarId(NbField)%size 
     1207            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),  & 
     1208                                  FieldVarId(NbField)%nc_id(i)) 
     1209            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat") 
     1210          ENDDO         
     1211        ENDIF  
     1212   
     1213       
     1214        status = NF90_ENDDEF(ncid)       
     1215 
     1216!        ncell=1 
     1217!        DO ind=ind_b,ind_e 
     1218!          d=>domain_type(ind) 
     1219   
     1220!          n=0 
     1221!          DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     1222!            DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     1223!              IF (single .OR. d%assign_domain(i,j)==ind) n=n+1 
     1224!            ENDDO 
     1225!          ENDDO 
     1226           
     1227         ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell)) 
     1228           
     1229         n=0   
     1230         DO ind=ind_b,ind_e 
     1231           d=>domain_type(ind) 
     1232           DO j=d%jj_begin-halo_size,d%jj_end+halo_size 
     1233             DO i=d%ii_begin-halo_size,d%ii_end+halo_size 
     1234               IF (d%assign_domain(i,j)==ind .OR. single) THEN 
     1235                 n=n+1 
     1236                 CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n)) 
     1237                 lon(n)=lon(n)*180/Pi 
     1238                 lat(n)=lat(n)*180/Pi 
     1239                 DO k=0,5 
     1240                   CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n)) 
     1241                   bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     1242                   bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     1243                 ENDDO 
     1244               ENDIF 
     1245             ENDDO 
     1246           ENDDO 
     1247         ENDDO 
     1248         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /)) 
     1249         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /)) 
     1250         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
     1251         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
     1252   
     1253!         ncell=ncell+n 
     1254         DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
     1255!      END DO  
     1256 
     1257    ELSE IF (Field(ind_b)%field_type==field_Z) THEN 
     1258        nvert=3 
     1259        DO ind=ind_b,ind_e 
     1260          d=>domain_type(ind) 
     1261         
     1262          DO j=d%jj_begin+1,d%jj_end 
     1263            DO i=d%ii_begin,d%ii_end-1 
     1264              ncell=ncell+1 
     1265            ENDDO 
     1266          ENDDO 
     1267 
     1268          DO j=d%jj_begin,d%jj_end-1 
     1269            DO i=d%ii_begin+1,d%ii_end 
     1270              ncell=ncell+1 
     1271            ENDDO 
     1272          ENDDO 
     1273 
     1274        END DO 
     1275       
     1276        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     1277        FieldId(NbField)=ncid 
     1278        status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     1279        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
     1280 
     1281        IF (Field(ind_b)%ndim==2)  THEN 
     1282          FieldVarId(NbField)%size=1 
     1283          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     1284        ELSE IF (Field(ind_b)%ndim==3)  THEN 
     1285          FieldVarId(NbField)%size=1 
     1286          ALLOCATE(FieldVarId(NbField)%nc_id(1)) 
     1287          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id) 
     1288        ELSE IF (Field(1)%ndim==4) THEN 
     1289          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3) 
     1290          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size)) 
     1291          status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id) 
     1292        ENDIF 
     1293 
     1294 
     1295       
     1296        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId) 
     1297       
     1298        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId) 
     1299        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude") 
     1300        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east") 
     1301        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon") 
     1302        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId) 
     1303        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude") 
     1304        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north") 
     1305        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat") 
     1306        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId) 
     1307        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId) 
     1308 
     1309 
     1310        IF (Field(ind_b)%ndim==2) THEN 
     1311          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
     1312          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1313        ELSE IF (Field(ind_b)%ndim==3) THEN 
     1314          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
     1315          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1316        ELSE IF (Field(ind_b)%ndim==4) THEN 
     1317          DO q=1,FieldVarId(NbField)%size 
     1318            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE,             & 
     1319                                  (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) 
     1320            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat") 
     1321          ENDDO         
     1322        ENDIF  
     1323         
     1324        status = NF90_ENDDEF(ncid)       
     1325 
     1326!        ncell=1 
     1327!        DO ind=ind_b,ind_e 
     1328!          d=>domain_type(ind) 
     1329!   
     1330!          n=0 
     1331!          DO j=d%jj_begin+1,d%jj_end 
     1332!            DO i=d%ii_begin,d%ii_end-1 
     1333!              n=n+1 
     1334!            ENDDO 
     1335!          ENDDO 
     1336! 
     1337!          DO j=d%jj_begin,d%jj_end-1 
     1338!            DO i=d%ii_begin+1,d%ii_end 
     1339!              n=n+1 
     1340!            ENDDO 
     1341!          ENDDO 
     1342 
     1343!         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n)) 
     1344         ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell)) 
     1345           
     1346         n=0   
     1347        
     1348         DO ind=ind_b,ind_e 
     1349           d=>domain_type(ind) 
     1350           DO j=d%jj_begin+1,d%jj_end 
     1351             DO i=d%ii_begin,d%ii_end-1 
     1352               nij=(j-1)*d%iim+i 
     1353               n=n+1 
     1354               CALL xyz2lonlat(d%vertex(:,vdown,i,j),lon(n),lat(n)) 
     1355               lon(n)=lon(n)*180/Pi 
     1356 !              IF (lon(n)<0) lon(n)=lon(n)+360 
     1357               lat(n)=lat(n)*180/Pi 
     1358!               CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
     1359!               CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
     1360!               CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
     1361 
     1362               CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n)) 
     1363               CALL xyz2lonlat(d%xyz(:,i,j-1),bounds_lon(1,n), bounds_lat(1,n)) 
     1364               CALL xyz2lonlat(d%xyz(:,i+1,j-1),bounds_lon(2,n), bounds_lat(2,n)) 
     1365 
     1366               DO k=0,2 
     1367                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     1368                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     1369!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
     1370               ENDDO 
     1371             ENDDO 
     1372           ENDDO 
     1373 
     1374           DO j=d%jj_begin,d%jj_end-1 
     1375             DO i=d%ii_begin+1,d%ii_end 
     1376               nij=(j-1)*d%iim+i 
     1377               n=n+1 
     1378!              CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n)) 
     1379               CALL xyz2lonlat(d%vertex(:,vup,i,j),lon(n),lat(n)) 
     1380               lon(n)=lon(n)*180/Pi 
     1381!              IF (lon(n)<0) lon(n)=lon(n)+360 
     1382               lat(n)=lat(n)*180/Pi 
     1383!              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
     1384!              CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n)) 
     1385!              CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n)) 
     1386               CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n)) 
     1387               CALL xyz2lonlat(d%xyz(:,i,j+1),bounds_lon(1,n), bounds_lat(1,n)) 
     1388               CALL xyz2lonlat(d%xyz(:,i-1,j+1),bounds_lon(2,n), bounds_lat(2,n)) 
     1389 
     1390               DO k=0,2 
     1391                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
     1392                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     1393!                 IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
     1394               ENDDO 
     1395             ENDDO 
     1396           ENDDO 
     1397         ENDDO  
     1398           
     1399         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /)) 
     1400         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /)) 
     1401         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
     1402         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /)) 
     1403   
     1404!          ncell=ncell+n 
     1405          DEALLOCATE(lon,lat,bounds_lon,bounds_lat) 
     1406!      END DO           
     1407    ENDIF 
     1408 
     1409 
     1410       
     1411!      status = NF90_CLOSE(ncid) 
     1412 
     1413    END SUBROUTINE Create_Header_gen 
     1414 
     1415    SUBROUTINE Create_header_mpi(name,field,nind) 
     1416    USE netcdf_mod 
    3041417    USE field_mod 
    3051418    USE domain_mod 
    … …  
    3071420    USE dimensions 
    3081421    USE geometry 
     1422    USE mpi_mod 
     1423    USE mpipara 
    3091424    IMPLICIT NONE 
    3101425      CHARACTER(LEN=*) :: name 
    … …  
    3241439      LOGICAL :: single  
    3251440      INTEGER :: nij 
     1441      INTEGER :: ncell_glo(0:mpi_size-1) 
     1442      INTEGER :: displ, ncell_tot 
     1443       
    3261444           
    3271445      NbField=NbField+1 
    … …  
    3571475        END DO 
    3581476       
    359         status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     1477        CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) 
     1478 
     1479        displ=0 
     1480        DO i=1,mpi_rank 
     1481          displ=displ+ncell_glo(i-1) 
     1482        ENDDO 
     1483        FieldVarId(NbField)%displ=displ 
     1484        ncell_tot=sum(ncell_glo(:)) 
     1485         
     1486        status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc', IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid) 
    3601487        FieldId(NbField)=ncid 
    361         status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     1488        status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId) 
    3621489        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
    3631490 
    … …  
    3921519          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
    3931520          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1521          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/)) 
    3941522        ELSE IF (Field(ind_b)%ndim==3) THEN 
    3951523          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
    3961524          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1525          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,size(field(ind_b)%rval3d,2),1/)) 
    3971526        ELSE IF (Field(ind_b)%ndim==4) THEN 
    3981527          DO i=1,FieldVarId(NbField)%size 
    399             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(i)) 
     1528            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),  & 
     1529                                  FieldVarId(NbField)%nc_id(i)) 
    4001530            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat") 
     1531            status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, (/ncell_tot,size(field(ind_b)%rval4d,2),1/)) 
    4011532          ENDDO         
    4021533        ENDIF  
    … …  
    4341565            ENDDO 
    4351566          ENDDO 
    436           status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) 
    437           status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) 
    438           status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
    439           status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     1567          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /)) 
     1568          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /)) 
     1569          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
     1570          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
    4401571   
    4411572          ncell=ncell+n 
    … …  
    4611592 
    4621593        END DO 
    463        
    464         status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
     1594         
     1595        CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr) 
     1596 
     1597        displ=0 
     1598        DO i=1,mpi_rank 
     1599          displ=displ+ncell_glo(i-1) 
     1600        ENDDO 
     1601        FieldVarId(NbField)%displ=displ 
     1602        ncell_tot=sum(ncell_glo(:)) 
     1603               
     1604        status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc',IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid) 
     1605!        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid) 
    4651606        FieldId(NbField)=ncid 
    466         status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId) 
     1607        status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId) 
    4671608        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid) 
    4681609 
    … …  
    4991640          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1)) 
    5001641          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1642          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/)) 
    5011643        ELSE IF (Field(ind_b)%ndim==3) THEN 
    5021644          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1)) 
    5031645          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat") 
     1646          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,size(field(ind_b)%rval3d,2),1/)) 
    5041647        ELSE IF (Field(ind_b)%ndim==4) THEN 
    5051648          DO q=1,FieldVarId(NbField)%size 
    506             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) 
     1649            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE,             & 
     1650                                  (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q)) 
    5071651            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat") 
     1652           status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, (/ncell_tot,size(field(ind_b)%rval4d,2),1/)) 
    5081653          ENDDO         
    5091654        ENDIF  
    … …  
    5751720           
    5761721           
    577           status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /)) 
    578           status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /)) 
    579           status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
    580           status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /)) 
     1722          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /)) 
     1723          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /)) 
     1724          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
     1725          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /)) 
    5811726   
    5821727          ncell=ncell+n 
    … …  
    5891734!      status = NF90_CLOSE(ncid) 
    5901735 
    591     end subroutine Create_Header 
     1736    end subroutine Create_Header_mpi   
    5921737    
    593    
     1738   SUBROUTINE Close_files 
     1739   USE netcdf 
     1740   IMPLICIT NONE 
     1741     INTEGER :: i,k,status 
     1742      
     1743     DO i=1,NbField 
     1744         status=NF90_CLOSE(FieldId(i)) 
     1745     ENDDO 
     1746   END SUBROUTINE  Close_files 
     1747      
     1748      
    5941749  function int2str(int) 
    5951750    implicit none 
Note: See TracChangeset for help on using the changeset viewer.