Changeset 15


Ignore:
Timestamp:
07/09/12 15:23:38 (12 years ago)
Author:
ymipsl
Message:

Update on 3D dynamic

YM

Location:
codes/icosagcm/trunk/src
Files:
3 added
1 deleted
16 edited

Legend:

Unmodified
Added
Removed

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

    r12 r15  
    1111  REAL(rstd),POINTER :: out_z(:,:) 
    1212  
     13  INTEGER :: itau_out 
     14   
    1315CONTAINS 
    14  
     16   
     17  SUBROUTINE init_caldyn(dt) 
     18  USE IOIPSL 
     19  IMPLICIT NONE 
     20    REAL(rstd),INTENT(IN) :: dt 
     21    REAL(rstd) :: write_period 
     22    CALL allocate_caldyn 
     23 
     24    CALL getin('write_period',write_period) 
     25     
     26    itau_out=INT(write_period/dt) 
     27   
     28  END SUBROUTINE init_caldyn 
     29  
    1530  SUBROUTINE allocate_caldyn 
    1631  USE domain_mod 
    … …  
    7893   
    7994   
    80   SUBROUTINE caldyn(f_phis, f_ps, f_theta_rhodz, f_u, f_dps, f_dtheta_rhodz, f_du) 
     95  SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_dps, f_dtheta_rhodz, f_du) 
    8196  USE domain_mod 
    8297  USE dimensions 
    … …  
    87102  USE vorticity_mod 
    88103  USE kinetic_mod 
     104  USE theta2theta_rhodz_mod 
    89105  IMPLICIT NONE 
     106  INTEGER,INTENT(IN)    :: it 
    90107  TYPE(t_field),POINTER :: f_phis(:) 
    91108  TYPE(t_field),POINTER :: f_ps(:) 
    … …  
    104121  REAL(rstd),POINTER :: du(:,:) 
    105122  INTEGER :: ind 
    106   INTEGER,SAVE :: it=0 
     123 
    107124   
    108125    CALL transfert_request(f_phis,req_i1)  
    … …  
    139156!    CALL kinetic(f_du,f_out) 
    140157 
    141     IF (mod(it,72)==0 ) THEN 
     158    IF (mod(it,itau_out)==0 ) THEN 
    142159      CALL writefield("ps",f_ps) 
    143       CALL writefield("dps",f_dps) 
     160!      CALL writefield("dps",f_dps) 
    144161!      CALL writefield("theta_rhodz",f_theta_rhodz) 
    145       CALL writefield("dtheta_rhodz",f_dtheta_rhodz) 
    146 !      CALL writefield("vort",f_out_z) 
    147       CALL writefield("theta",f_out) 
     162!      CALL writefield("dtheta_rhodz",f_dtheta_rhodz) 
     163      CALL vorticity(f_u,f_out_z) 
     164      CALL writefield("vort",f_out_z) 
     165!      CALL writefield("theta",f_out) 
     166      CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_out) ; 
     167      CALL writefield("T",f_out) 
    148168     
    149169!      CALL writefield("out",f_out) 
    … …  
    155175    ENDIF 
    156176!    CALL check_mass_conservation(f_ps,f_dps) 
    157    it=it+1       
     177 
    158178  END SUBROUTINE caldyn 
    159179   

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

    r12 r15  
    170170      CALL writefield("dh",f_dh) 
    171171      CALL Compute_enstrophy 
     172!      CALL Compute_PV 
    172173    ENDIF 
    173174    it=it+1       

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

    r12 r15  
    190190            ng2=d%neighbour(:,MOD(k+1,6),i,j) 
    191191            IF (sqrt(sum((ng1-ng2)**2))<1e-16) ng2=d%neighbour(:,MOD(k+2,6),i,j) 
    192             CALL circonscrit(d%xyz(:,i,j),ng1,ng2,d%vertex(:,k,i,j)) 
     192            CALL circumcenter(d%xyz(:,i,j),ng1,ng2,d%vertex(:,k,i,j)) 
    193193          ENDDO 
    194194        ENDDO 

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

    r12 r15  
    140140  REAL(rstd) :: r2 
    141141  REAL(rstd) :: utot 
     142  REAL(rstd) :: lonx,latx 
    142143  
    143144    DO l=1,llm 
    … …  
    157158     
    158159    CALL lonlat2xyz(Pi/9,2*Pi/9,V0) 
     160 
    159161    u(:,:)=1e10       
    160162    DO l=1,llm 
    … …  
    165167          CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat) 
    166168          CALL cross_product2(V0,xyz_e(ij+u_right,:)/radius,ep) 
    167           r2=(asin(sum(ep*ep)))**2 
     169          r2=(asin(sqrt(sum(ep*ep))))**2 
    168170          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 
    169171 
    … …  
    184186          CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat) 
    185187          CALL cross_product2(V0,xyz_e(ij+u_lup,:)/radius,ep) 
    186           r2=(asin(sum(ep*ep)))**2 
     188          r2=(asin(sqrt(sum(ep*ep))))**2 
    187189          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 
    188190          ulon(1) = -sin(lon) * utot 
    … …  
    205207          CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat) 
    206208          CALL cross_product2(V0,xyz_e(ij+u_ldown,:)/radius,ep) 
    207           r2=(asin(sum(ep*ep)))**2 
     209          r2=(asin(sqrt(sum(ep*ep))))**2 
    208210          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-r2/0.01) 
    209211          ulon(1) = -sin(lon) * utot 

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

    r12 r15  
    44  TYPE t_geometry 
    55    TYPE(t_field),POINTER :: xyz_i(:) 
     6    TYPE(t_field),POINTER :: centroid(:) 
    67    TYPE(t_field),POINTER :: xyz_e(:) 
    78    TYPE(t_field),POINTER :: xyz_v(:) 
     9    TYPE(t_field),POINTER :: ep_e(:) 
     10    TYPE(t_field),POINTER :: et_e(:) 
     11    TYPE(t_field),POINTER :: elon_i(:) 
     12    TYPE(t_field),POINTER :: elat_i(:) 
     13    TYPE(t_field),POINTER :: elon_e(:) 
     14    TYPE(t_field),POINTER :: elat_e(:) 
    815    TYPE(t_field),POINTER :: Ai(:) 
    916    TYPE(t_field),POINTER :: Av(:) 
    … …  
    2027  TYPE(t_geometry),TARGET :: geom 
    2128   
    22   REAL(rstd),POINTER :: Ai(:) 
    23   REAL(rstd),POINTER :: xyz_i(:,:) 
    24   REAL(rstd),POINTER :: xyz_e(:,:) 
    25   REAL(rstd),POINTER :: xyz_v(:,:) 
    26   REAL(rstd),POINTER :: Av(:) 
    27   REAL(rstd),POINTER :: de(:) 
    28   REAL(rstd),POINTER :: le(:) 
    29   REAL(rstd),POINTER :: Riv(:,:) 
    30   REAL(rstd),POINTER :: Riv2(:,:) 
    31   INTEGER,POINTER :: ne(:,:) 
    32   REAL(rstd),POINTER :: Wee(:,:,:) 
    33   REAL(rstd),POINTER :: bi(:) 
    34   REAL(rstd),POINTER :: fv(:) 
     29  REAL(rstd),POINTER :: Ai(:)          ! area of a cell 
     30  REAL(rstd),POINTER :: xyz_i(:,:)     ! coordinate of the center of the cell (voronoi) 
     31  REAL(rstd),POINTER :: centroid(:,:)  ! coordinate of the centroid of the cell 
     32  REAL(rstd),POINTER :: xyz_e(:,:)     ! coordinate of a wind point on the cell on a edge 
     33  REAL(rstd),POINTER :: ep_e(:,:)      ! perpendicular unit vector of a edge (outsider) 
     34  REAL(rstd),POINTER :: et_e(:,:)      ! tangeantial unit vector of a edge 
     35  REAL(rstd),POINTER :: elon_i(:,:)    ! unit longitude vector on the center  
     36  REAL(rstd),POINTER :: elat_i(:,:)    ! unit latitude vector on the center  
     37  REAL(rstd),POINTER :: elon_e(:,:)    ! unit longitude vector on a wind point  
     38  REAL(rstd),POINTER :: elat_e(:,:)    ! unit latitude vector on a wind point 
     39  REAL(rstd),POINTER :: xyz_v(:,:)     ! coordinate of a vertex (center of the dual mesh) 
     40  REAL(rstd),POINTER :: Av(:)          ! area of dual mesk cell 
     41  REAL(rstd),POINTER :: de(:)          ! distance from a neighbour == lenght of an edge of the dual mesh 
     42  REAL(rstd),POINTER :: le(:)          ! lenght of a edge 
     43  REAL(rstd),POINTER :: Riv(:,:)       ! weight 
     44  REAL(rstd),POINTER :: Riv2(:,:)      ! weight 
     45  INTEGER,POINTER    :: ne(:,:)        ! convention for the way on the normal wind on an edge  
     46  REAL(rstd),POINTER :: Wee(:,:,:)     ! weight 
     47  REAL(rstd),POINTER :: bi(:)          ! orographie 
     48  REAL(rstd),POINTER :: fv(:)          ! coriolis (evaluted on a vertex) 
    3549      
    3650 
    … …  
    4458    CALL allocate_field(geom%Ai,field_t,type_real) 
    4559    CALL allocate_field(geom%xyz_i,field_t,type_real,3) 
     60    CALL allocate_field(geom%centroid,field_t,type_real,3) 
    4661    CALL allocate_field(geom%xyz_e,field_u,type_real,3) 
     62    CALL allocate_field(geom%ep_e,field_u,type_real,3) 
     63    CALL allocate_field(geom%et_e,field_u,type_real,3) 
     64    CALL allocate_field(geom%elon_i,field_t,type_real,3) 
     65    CALL allocate_field(geom%elat_i,field_t,type_real,3) 
     66    CALL allocate_field(geom%elon_e,field_u,type_real,3) 
     67    CALL allocate_field(geom%elat_e,field_u,type_real,3) 
    4768    CALL allocate_field(geom%xyz_v,field_z,type_real,3) 
    4869    CALL allocate_field(geom%de,field_u,type_real) 
    … …  
    6687    Ai=geom%Ai(ind) 
    6788    xyz_i=geom%xyz_i(ind) 
     89    centroid=geom%centroid(ind) 
    6890    xyz_e=geom%xyz_e(ind) 
     91    ep_e=geom%ep_e(ind) 
     92    et_e=geom%et_e(ind) 
     93    elon_i=geom%elon_i(ind) 
     94    elat_i=geom%elat_i(ind) 
     95    elon_e=geom%elon_e(ind) 
     96    elat_e=geom%elat_e(ind) 
    6997    xyz_v=geom%xyz_v(ind) 
    7098    de=geom%de(ind) 
    … …  
    79107     
    80108  END SUBROUTINE swap_geometry 
    81      
    82   SUBROUTINE set_geometry 
     109   
     110  SUBROUTINE optimize_geometry 
    83111  USE metric 
    84112  USE spherical_geom_mod 
    85113  USE domain_mod 
    86114  USE dimensions 
     115  USE transfert_mod 
     116  USE vector 
     117  IMPLICIT NONE 
     118    INTEGER,PARAMETER :: nb_it=3000 
     119    TYPE(t_domain),POINTER :: d 
     120    INTEGER :: ind,it,i,j,n,k 
     121    REAL(rstd) :: x1(3),x2(3) 
     122    REAL(rstd) :: vect(3,6) 
     123    REAL(rstd) :: centr(3) 
     124    REAL(rstd) :: sum 
     125    LOGICAL    :: check 
     126     
     127    DO ind=1,ndomain 
     128      d=>domain(ind) 
     129      CALL swap_dimensions(ind) 
     130      CALL swap_geometry(ind) 
     131      DO j=jj_begin,jj_end 
     132        DO i=ii_begin,ii_end 
     133          n=(j-1)*iim+i 
     134          xyz_i(n,:)=d%xyz(:,i,j)  
     135        ENDDO 
     136      ENDDO 
     137    ENDDO 
     138     
     139    CALL transfert_request(geom%xyz_i,req_i1) 
     140     
     141    DO ind=1,ndomain 
     142      CALL swap_dimensions(ind) 
     143      CALL swap_geometry(ind) 
     144      DO j=jj_begin,jj_end 
     145        DO i=ii_begin,ii_end 
     146          n=(j-1)*iim+i 
     147          DO k=0,5 
     148            x1(:) = xyz_i(n+t_pos(k+1),:) 
     149            x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:) 
     150            if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:) 
     151            CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:)) 
     152          ENDDO 
     153        ENDDO 
     154      ENDDO 
     155    ENDDO 
     156     
     157    DO ind=1,ndomain 
     158      d=>domain(ind) 
     159      CALL swap_dimensions(ind) 
     160      CALL swap_geometry(ind) 
     161      DO j=jj_begin,jj_end 
     162        DO i=ii_begin,ii_end 
     163          n=(j-1)*iim+i 
     164          DO k=0,5 
     165            x1(:) = xyz_v(n+z_pos(k+1),:) 
     166            x2(:) = d%vertex(:,k,i,j)  
     167            IF (norm(x1-x2)>1e-10) THEN 
     168              PRINT*,"vertex diff ",ind,i,j,k 
     169              PRINT*,x1 
     170              PRINT*,x2 
     171            ENDIF 
     172          ENDDO 
     173        ENDDO 
     174      ENDDO 
     175    ENDDO 
     176     
     177     
     178    DO it=1,nb_it 
     179      IF (MOD(it,100)==0) THEN  
     180        check=.TRUE. 
     181      ELSE 
     182        check=.FALSE. 
     183      ENDIF 
     184       
     185      sum=0 
     186      DO ind=1,ndomain 
     187        CALL swap_dimensions(ind) 
     188        CALL swap_geometry(ind) 
     189        DO j=jj_begin,jj_end 
     190          DO i=ii_begin,ii_end 
     191            n=(j-1)*iim+i 
     192            vect(:,1)=xyz_v(n+z_rup,:) 
     193            vect(:,2)=xyz_v(n+z_up,:) 
     194            vect(:,3)=xyz_v(n+z_lup,:) 
     195            vect(:,4)=xyz_v(n+z_ldown,:) 
     196            vect(:,5)=xyz_v(n+z_down,:) 
     197            vect(:,6)=xyz_v(n+z_rdown,:) 
     198            CALL compute_centroid(vect,6,centr) 
     199            IF (check) THEN 
     200              sum=MAX(sum,norm(xyz_i(n,:)-centr(:))) 
     201            ENDIF 
     202            xyz_i(n,:)=centr(:) 
     203          ENDDO 
     204        ENDDO 
     205!        i=ii_begin ; j=jj_begin ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j)  
     206!        i=ii_begin ; j=jj_end   ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j) 
     207!        i=ii_end   ; j=jj_begin ; n=(j-1)*iim+i ;   xyz_i(n,:)=domain(ind)%xyz(:,i,j) 
     208!        PRINT *,"Pb ?? : "  
     209!        PRINT *, xyz_i(n,:), domain(ind)%xyz(:,i,j), norm(xyz_i(n,:)- domain(ind)%xyz(:,i,j)) 
     210!        i=ii_end   ; j=jj_end   ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j) 
     211         
     212      ENDDO 
     213 
     214      IF (check) THEN 
     215!        sum=sum/(ndomain*ii_nb*jj_nb) 
     216        PRINT *,"it = ",it,"  diff centroid circumcenter ",sum 
     217      ENDIF  
     218       
     219     CALL transfert_request(geom%xyz_i,req_i1) 
     220 
     221    DO ind=1,ndomain 
     222      CALL swap_dimensions(ind) 
     223      CALL swap_geometry(ind) 
     224      DO j=jj_begin,jj_end 
     225        DO i=ii_begin,ii_end 
     226          n=(j-1)*iim+i 
     227          DO k=0,5 
     228            x1(:) = xyz_i(n+t_pos(k+1),:) 
     229            x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:) 
     230            if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:) 
     231            CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:)) 
     232          ENDDO 
     233        ENDDO 
     234      ENDDO 
     235    ENDDO 
     236 
     237       
     238    ENDDO 
     239     
     240  END SUBROUTINE optimize_geometry 
     241   
     242     
     243  SUBROUTINE set_geometry 
     244  USE metric 
     245  USE vector 
     246  USE spherical_geom_mod 
     247  USE domain_mod 
     248  USE dimensions 
     249  USE transfert_mod 
    87250  IMPLICIT NONE 
    88251 
    89252    REAL(rstd) :: surf(6) 
    90253    REAL(rstd) :: surf_v(6) 
     254    REAL(rstd) :: vect(3,6) 
     255    REAL(rstd) :: centr(6) 
     256    REAL(rstd) :: vet(3),vep(3) 
    91257    INTEGER :: ind,i,j,k,n 
    92258    TYPE(t_domain),POINTER :: d 
    93259    REAL(rstd) ::  S 
    94260    REAL(rstd) :: w(6) 
    95     REAl(rstd) :: lon,lat 
     261    REAL(rstd) :: lon,lat 
    96262    INTEGER :: ii_glo,jj_glo 
    97263    REAL(rstd) :: S1,S2 
     264           
    98265       
     266    CALL optimize_geometry 
     267     
    99268    DO ind=1,ndomain 
    100269      d=>domain(ind) 
    … …  
    104273        DO i=ii_begin-1,ii_end+1 
    105274          n=(j-1)*iim+i 
    106          
    107           xyz_i(n,:)=d%xyz(:,i,j)  
    108           xyz_v(n+z_up,:)=d%vertex(:,vup-1,i,j)  
    109           xyz_v(n+z_down,:)=d%vertex(:,vdown-1,i,j)  
    110            
     275 
     276          DO k=0,5 
     277            ne(n,k+1)=d%ne(k,i,j) 
     278          ENDDO 
     279                   
     280!          xyz_i(n,:)=d%xyz(:,i,j)  
     281!          xyz_v(n+z_up,:)=d%vertex(:,vup-1,i,j)  
     282!          xyz_v(n+z_down,:)=d%vertex(:,vdown-1,i,j)  
     283 
     284          vect(:,1)=xyz_v(n+z_rup,:) 
     285          vect(:,2)=xyz_v(n+z_up,:) 
     286          vect(:,3)=xyz_v(n+z_lup,:) 
     287          vect(:,4)=xyz_v(n+z_ldown,:) 
     288          vect(:,5)=xyz_v(n+z_down,:) 
     289          vect(:,6)=xyz_v(n+z_rdown,:) 
     290          CALL compute_centroid(vect,6,centr) 
     291          centroid(n,:)=centr(:) 
     292 
     293                
    111294          CALL xyz2lonlat(xyz_v(n+z_up,:),lon,lat) 
    112295          fv(n+z_up)=2*sin(lat)*omega 
    … …  
    116299          bi(n)=0.  
    117300         
    118           CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),de(n+u_right)) 
    119           CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),de(n+u_lup)) 
    120           CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),de(n+u_ldown)) 
    121            
    122           CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),0.5,xyz_e(n+u_right,:)) 
    123           CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),0.5,xyz_e(n+u_lup,:)) 
    124           CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),0.5,xyz_e(n+u_ldown,:)) 
    125            
    126           CALL dist_cart(d%vertex(:,vrdown-1,i,j),d%vertex(:,vrup-1,i,j),le(n+u_right)) 
    127           CALL dist_cart(d%vertex(:,vup-1,i,j),d%vertex(:,vlup-1,i,j),le(n+u_lup)) 
    128           CALL dist_cart(d%vertex(:,vldown-1,i,j),d%vertex(:,vdown-1,i,j),le(n+u_ldown)) 
    129  
     301!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),de(n+u_right)) 
     302!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),de(n+u_lup)) 
     303!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),de(n+u_ldown)) 
     304 
     305          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_right,:),de(n+u_right)) 
     306          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_lup,:),de(n+u_lup)) 
     307          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_ldown,:),de(n+u_ldown)) 
     308           
     309!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),0.5,xyz_e(n+u_right,:)) 
     310!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),0.5,xyz_e(n+u_lup,:)) 
     311!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),0.5,xyz_e(n+u_ldown,:)) 
     312 
     313          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_right,:),0.5,xyz_e(n+u_right,:)) 
     314          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_lup,:),0.5,xyz_e(n+u_lup,:)) 
     315          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_ldown,:),0.5,xyz_e(n+u_ldown,:)) 
     316          
     317!          CALL dist_cart(d%vertex(:,vrdown-1,i,j),d%vertex(:,vrup-1,i,j),le(n+u_right)) 
     318!          CALL dist_cart(d%vertex(:,vup-1,i,j),d%vertex(:,vlup-1,i,j),le(n+u_lup)) 
     319!          CALL dist_cart(d%vertex(:,vldown-1,i,j),d%vertex(:,vdown-1,i,j),le(n+u_ldown)) 
     320 
     321          CALL dist_cart(xyz_v(n+z_rdown,:), xyz_v(n+z_rup,:),le(n+u_right)) 
     322          CALL dist_cart(xyz_v(n+z_up,:), xyz_v(n+z_lup,:),le(n+u_lup)) 
     323          CALL dist_cart(xyz_v(n+z_ldown,:), xyz_v(n+z_down,:),le(n+u_ldown)) 
    130324          
    131325          Ai(n)=0 
    132326          DO k=0,5 
    133             CALL surf_triangle(d%xyz(:,i,j),d%neighbour(:,k,i,j),d%neighbour(:,MOD((k+1+6),6),i,j),surf_v(k+1)) 
    134             CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,MOD((k-1+6),6),i,j),d%vertex(:,k,i,j),surf(k+1)) 
    135             ne(n,k+1)=d%ne(k,i,j) 
     327!            CALL surf_triangle(d%xyz(:,i,j),d%neighbour(:,k,i,j),d%neighbour(:,MOD((k+1+6),6),i,j),surf_v(k+1)) 
     328!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,MOD((k-1+6),6),i,j),d%vertex(:,k,i,j),surf(k+1)) 
     329            CALL surf_triangle(xyz_i(n,:),xyz_i(n+t_pos(k+1),:),xyz_i(n+t_pos(MOD((k+1+6),6)+1),:),surf_v(k+1)) 
     330            CALL surf_triangle(xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:),surf(k+1)) 
    136331            Ai(n)=Ai(n)+surf(k+1) 
    137           ENDDO 
    138  
    139  
     332            IF (i==ii_end .AND. j==jj_begin) THEN 
     333              IF (Ai(n)<1e20) THEN 
     334              ELSE 
     335                PRINT *,"PB !!",Ai(n),k,surf(k+1) 
     336                PRINT*,xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:) 
     337              ENDIF 
     338            ENDIF 
     339          ENDDO 
     340 
     341 
     342          CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),vep) 
     343          IF (norm(vep)>1e-30) THEN 
     344            vep(:)=vep(:)/norm(vep) 
     345            CALL cross_product2(vep,xyz_e(n+u_right,:),vet) 
     346            vet(:)=vet(:)/norm(vet) 
     347            ep_e(n+u_right,:)=vep(:)*ne(n,right) 
     348            et_e(n+u_right,:)=vet(:)*ne(n,right) 
     349          ENDIF 
     350 
     351          CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),vep) 
     352          IF (norm(vep)>1e-30) THEN 
     353            vep(:)=vep(:)/norm(vep) 
     354            CALL cross_product2(vep,xyz_e(n+u_lup,:),vet) 
     355            vet(:)=vet(:)/norm(vet) 
     356            ep_e(n+u_lup,:)=vep(:)*ne(n,lup) 
     357            et_e(n+u_lup,:)=vet(:)*ne(n,lup) 
     358          ENDIF 
     359 
     360          CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),vep) 
     361          IF (norm(vep)>1e-30) THEN 
     362            vep(:)=vep(:)/norm(vep) 
     363            CALL cross_product2(vep,xyz_e(n+u_ldown,:),vet) 
     364            vet(:)=vet(:)/norm(vet) 
     365            ep_e(n+u_ldown,:)=vep(:)*ne(n,ldown) 
     366            et_e(n+u_ldown,:)=vet(:)*ne(n,ldown) 
     367          ENDIF 
     368 
     369          CALL xyz2lonlat(xyz_i(n,:),lon,lat) 
     370          elon_i(n,1) = -sin(lon) 
     371          elon_i(n,2) = cos(lon) 
     372          elon_i(n,3) = 0 
     373          elat_i(n,1) = -cos(lon)*sin(lat) 
     374          elat_i(n,2) = -sin(lon)*sin(lat) 
     375          elat_i(n,3) = cos(lat) 
     376 
     377           
     378          CALL xyz2lonlat(xyz_e(n+u_right,:),lon,lat) 
     379          elon_e(n+u_right,1) = -sin(lon) 
     380          elon_e(n+u_right,2) = cos(lon) 
     381          elon_e(n+u_right,3) = 0 
     382          elat_e(n+u_right,1) = -cos(lon)*sin(lat) 
     383          elat_e(n+u_right,2) = -sin(lon)*sin(lat) 
     384          elat_e(n+u_right,3) = cos(lat) 
     385           
     386          CALL xyz2lonlat(xyz_e(n+u_lup,:),lon,lat) 
     387          elon_e(n+u_lup,1) = -sin(lon) 
     388          elon_e(n+u_lup,2) = cos(lon) 
     389          elon_e(n+u_lup,3) = 0 
     390          elat_e(n+u_lup,1) = -cos(lon)*sin(lat) 
     391          elat_e(n+u_lup,2) = -sin(lon)*sin(lat) 
     392          elat_e(n+u_lup,3) = cos(lat) 
     393  
     394          CALL xyz2lonlat(xyz_e(n+u_ldown,:),lon,lat) 
     395          elon_e(n+u_ldown,1) = -sin(lon) 
     396          elon_e(n+u_ldown,2) = cos(lon) 
     397          elon_e(n+u_ldown,3) = 0 
     398          elat_e(n+u_ldown,1) = -cos(lon)*sin(lat) 
     399          elat_e(n+u_ldown,2) = -sin(lon)*sin(lat) 
     400          elat_e(n+u_ldown,3) = cos(lat) 
     401  
     402           
    140403          DO k=0,5 
    141             CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,k,i,j),S1) 
    142             CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,MOD(k+1+6,6),i,j),S2) 
     404!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,k,i,j),S1) 
     405!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,MOD(k+1+6,6),i,j),S2) 
     406            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1) 
     407            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(MOD(k+1+6,6)+1),:),S2) 
    143408            Riv(n,k+1)=0.5*(S1+S2)/Ai(n) 
    144409            Riv2(n,k+1)=0.5*(S1+S2)/surf_v(k+1) 
    … …  
    236501                  
    237502    ENDDO 
     503    CALL transfert_request(geom%Ai,req_i1) 
     504    CALL transfert_request(geom%centroid,req_i1) 
    238505    CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S) 
    239506!    PRINT *,"Surf triangle : ",S*20/(4*Pi) 

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

    r12 r15  
    11MODULE grid_param 
    2   INTEGER,PARAMETER  :: iim_glo=40 
    3   INTEGER,PARAMETER  :: jjm_glo=iim_glo 
     2  INTEGER  :: iim_glo=40 
     3  INTEGER  :: jjm_glo 
    44  INTEGER,PARAMETER  :: nb_face=10 
    5   INTEGER,PARAMETER  :: llm=19 
     5  INTEGER  :: llm=19 
     6 
     7CONTAINS 
     8 
     9  SUBROUTINE init_grid_param 
     10  USE IOIPSL 
     11  IMPLICIT NONE 
     12    CALL getin('nbp',iim_glo) 
     13    jjm_glo=40 
     14    CALL getin('llm',llm) 
     15     
     16  END SUBROUTINE  init_grid_param 
     17 
    618END MODULE grid_param 
    719     

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

    r12 r15  
    99  USE timeloop_gcm_mod 
    1010  USE transfert_mod 
    11   USE dissip_mod 
    1211  USE disvert_mod 
    1312  USE etat0_mod 
    1413  USE transfert_mod 
     14  USE vector 
     15  USE wind_mod 
     16  USE grid_param 
    1517  IMPLICIT NONE 
    1618   
    … …  
    2022  INTEGER :: ind,i,j,k,n 
    2123  REAL(rstd) :: tot_sum=0 
     24  REAL(rstd) :: vect(3,6) 
     25  REAL(rstd) :: centr(3),dist 
    2226  
     27  CALL init_grid_param 
    2328  CALL compute_metric 
    2429  CALL compute_domain 
     30  CALL init_transfert 
    2531  CALL compute_geometry 
    26   CALL init_transfert 
    2732  CALL init_disvert  
    2833   
    … …  
    5459 
    5560   
    56 !  CALL WriteField("Ai",geom%Ai) 
     61  CALL WriteField("Ai",geom%Ai) 
    5762!  CALL WriteField("sum_ne",sum_ne) 
    58  
    59  
     63   
    6064  CALL timeloop 
    6165   

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

    r13 r15  
    99  USE timeloop_sw_mod 
    1010  USE transfert_mod 
    11   USE dissip_mod 
     11  USE dissip_sw_mod 
    1212  USE disvert_mod 
    1313  USE transfert_mod 

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

    r13 r15  
    128128!  PRINT *,"dist",vertex_glo(2,1,1)%xyz 
    129129!  PRINT *,"dist",p1/sqrt(sum(p1**2)) 
    130   CALL circonscrit(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,2,1)%xyz,p1) 
     130  CALL circumcenter(vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,2,1)%xyz,p1) 
    131131!  CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(2,1,1)%xyz,vertex_glo(1,1,1)%xyz,p1) 
    132132!  CALL Centroide(vertex_glo(1,2,1)%xyz,vertex_glo(1,1,1)%xyz,vertex_glo(2,1,1)%xyz,p1) 

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

    r12 r15  
    3939 
    4040    DO    l    = 1, llm+1 
     41!$OMP DO 
    4142      DO j=jj_begin-offset,jj_end+offset 
    4243        DO i=ii_begin-offset,ii_end+offset 

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

    r13 r15  
    7878 
    7979  REAL(rstd)  :: AB,AC,BC 
    80   REAL(rstd)  :: s 
     80  REAL(rstd)  :: s,x 
    8181   
    8282  CALL dist_cart(A,B,AB) 
    … …  
    8585   
    8686  s=(AB+AC+BC)/2 
    87   surf=4*atan(sqrt( tan(s/2) * tan((s-AB)/2)  * tan((s-AC)/2) * tan((s-BC)/2))) 
     87  x=tan(s/2) * tan((s-AB)/2)  * tan((s-AC)/2) * tan((s-BC)/2) 
     88  IF (x<0) x=0. 
     89  surf=4*atan(sqrt( x)) 
    8890   
    8991END SUBROUTINE surf_triangle  
    … …  
    142144 
    143145 
    144   SUBROUTINE circonscrit(A0,B0,C0,Center) 
     146  SUBROUTINE circumcenter(A0,B0,C0,Center) 
    145147  USE vector 
    146148  IMPLICIT NONE 
    … …  
    157159    center=center/sqrt(sum(center**2)) 
    158160     
    159   END SUBROUTINE circonscrit 
     161  END SUBROUTINE circumcenter 
    160162     
    161163 
    162   SUBROUTINE centroide(A,B,C,Center) 
     164  SUBROUTINE compute_centroid(points,n,centr) 
     165  USE vector 
    163166  IMPLICIT NONE 
    164     REAL(rstd), INTENT(IN)  :: A(3),B(3),C(3) 
    165     REAL(rstd), INTENT(OUT) :: Center(3) 
     167    INTEGER :: n 
     168    REAL(rstd), INTENT(IN)  :: points(3,n) 
     169    REAL(rstd), INTENT(OUT) :: Centr(3) 
     170     
     171    REAL(rstd) :: p1(3),p2(3),cross(3) 
     172    REAL(rstd) :: norm_cross 
     173    INTEGER :: i,j 
     174       
     175      Centr(:)=0 
     176      DO i=1,n 
     177        j=MOD(i,n)+1 
     178        p1=points(:,i)/norm(points(:,i)) 
     179        p2=points(:,j)/norm(points(:,j)) 
     180        CALL cross_product2(p1,p2,cross) 
     181        norm_cross=norm(cross) 
     182        if (norm_cross<1e-10) CYCLE  
     183         
     184        Centr(:)=centr(:)+asin(norm_cross)*cross(:)/norm_cross 
     185      ENDDO 
     186       
     187      Centr(:)=centr(:)/norm(centr(:)) 
    166188   
    167     REAL(rstd)  :: AB05(3),AC05(3),BC05(3) 
    168     REAL(rstd)  :: Vab(3),Vbc(3),Vca(3) 
    169     REAL(rstd)  :: t 
    170     INTEGER :: n1,n2 
    171         
    172     AB05=(A+B)/2. 
    173     AC05=(A+C)/2. 
    174     BC05=(B+C)/2. 
    175      
    176     CALL director(A,B,C,Vab) 
    177 !    PRINT *,'director',sum((B-A)*Vab) 
    178     CALL director(B,C,A,Vbc) 
    179 !    PRINT *,'director',sum((C-B)*Vbc) 
    180     CALL director(C,A,B,Vca) 
    181 !    PRINT *,'director',sum((A-C)*Vca) 
    182     
    183     n1=1 ; n2=2 
    184      
    185     t=(Vbc(n1)*(BC05(n2)-AB05(n2))+(AB05(n1)-BC05(n1))*Vbc(n2))/(Vbc(n1)*Vab(n2)-Vab(n1)*Vbc(n2)) 
    186     Center=AB05+t*Vab 
    187     Center=Center/sqrt(sum(Center**2)) 
    188 !    PRINT*,"Center :",A 
    189 !    PRINT*,"Center :",B 
    190 !    PRINT*,"Center :",C 
    191 !    PRINT*,"Center :",Center 
    192     
    193   CONTAINS  
    194    
    195     SUBROUTINE director(A,B,C,V) 
    196     IMPLICIT NONE 
    197     REAL(rstd), INTENT(IN)  :: A(3),B(3),C(3) 
    198     REAL(rstd), INTENT(OUT) :: V(3) 
    199     REAL(rstd) :: AB(3),AC(3) 
    200     REAL(rstd) :: t 
    201      
    202       AB=B-A 
    203       AC=C-A 
    204      
    205       t=sum(AB*AC)/sum(AB**2) 
    206       V=AB*t-AC 
    207     END SUBROUTINE director 
    208         
    209   END SUBROUTINE centroide 
     189  END SUBROUTINE compute_centroid 
    210190 
    211191END MODULE spherical_geom_mod 

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

    r12 r15  
    22 
    33CONTAINS 
     4   
     5  SUBROUTINE theta_rhodz2theta(f_ps,f_theta_rhodz,f_theta) 
     6  USE transfert_mod 
     7  USE field_mod 
     8  USE dimensions 
     9  USE geometry 
     10  USE domain_mod 
     11  IMPLICIT NONE 
     12    TYPE(t_field), POINTER :: f_ps(:) 
     13    TYPE(t_field), POINTER :: f_theta_rhodz(:) 
     14    TYPE(t_field), POINTER :: f_theta(:) 
     15   
     16    REAL(rstd), POINTER :: ps(:) 
     17    REAL(rstd), POINTER :: theta_rhodz(:,:) 
     18    REAL(rstd), POINTER :: theta(:,:) 
     19    INTEGER :: ind 
    420 
     21    DO ind=1,ndomain 
     22      CALL swap_dimensions(ind) 
     23      CALL swap_geometry(ind) 
     24      ps=f_ps(ind) 
     25      theta_rhodz=f_theta_rhodz(ind) 
     26      theta=f_theta(ind) 
     27      CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0) 
     28    ENDDO 
     29   
     30  END SUBROUTINE theta_rhodz2theta 
     31 
     32  SUBROUTINE theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp) 
     33  USE transfert_mod 
     34  USE field_mod 
     35  USE dimensions 
     36  USE geometry 
     37  USE domain_mod 
     38  IMPLICIT NONE 
     39    TYPE(t_field), POINTER :: f_ps(:) 
     40    TYPE(t_field), POINTER :: f_theta_rhodz(:) 
     41    TYPE(t_field), POINTER :: f_temp(:) 
     42   
     43    REAL(rstd), POINTER :: ps(:) 
     44    REAL(rstd), POINTER :: theta_rhodz(:,:) 
     45    REAL(rstd), POINTER :: temp(:,:) 
     46    INTEGER :: ind 
     47 
     48    DO ind=1,ndomain 
     49      CALL swap_dimensions(ind) 
     50      CALL swap_geometry(ind) 
     51      ps=f_ps(ind) 
     52      theta_rhodz=f_theta_rhodz(ind) 
     53      temp=f_temp(ind) 
     54      CALL compute_theta_rhodz2temperature(ps, theta_rhodz,temp,0) 
     55    ENDDO 
     56   
     57  END SUBROUTINE theta_rhodz2temperature 
     58     
    559  SUBROUTINE theta2theta_rhodz(f_ps,f_theta,f_theta_rhodz) 
    660  USE transfert_mod 
    … …  
    4195    INTEGER,INTENT(IN) :: offset 
    4296    INTEGER :: i,j,ij,l 
    43     REAL(rstd) :: p(iim*jjm,llm+1) 
     97    REAL(rstd),ALLOCATABLE,SAVE :: p(:,:) 
     98 
     99!$OMP BARRIER 
     100!$OMP MASTER 
     101      ALLOCATE( p(iim*jjm,llm+1)) 
     102!$OMP END MASTER 
     103!$OMP BARRIER 
    44104 
    45105    CALL compute_pression(ps,p,offset) 
    46106     
    47107    DO    l    = 1, llm 
     108!$OMP DO 
    48109      DO j=jj_begin-offset,jj_end+offset 
    49110        DO i=ii_begin-offset,ii_end+offset 
    … …  
    53114      ENDDO 
    54115    ENDDO 
    55      
     116 
     117!$OMP BARRIER 
     118!$OMP MASTER 
     119      DEALLOCATE( p) 
     120!$OMP END MASTER 
     121!$OMP BARRIER     
     122 
    56123  END SUBROUTINE compute_theta2theta_rhodz 
    57124 
     125  SUBROUTINE compute_theta_rhodz2theta(ps,theta_rhodz,theta,offset) 
     126  USE dimensions 
     127  USE geometry 
     128  USE metric 
     129  USE pression_mod 
     130  IMPLICIT NONE 
     131    REAL(rstd),INTENT(IN) :: ps(iim*jjm) 
     132    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm) 
     133    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 
     134    INTEGER,INTENT(IN) :: offset 
     135    INTEGER :: i,j,ij,l 
     136    REAL(rstd),SAVE,ALLOCATABLE :: p(:,:) 
     137 
     138!$OMP BARRIER 
     139!$OMP MASTER 
     140      ALLOCATE( p(iim*jjm,llm+1)) 
     141!$OMP END MASTER 
     142!$OMP BARRIER 
     143    CALL compute_pression(ps,p,offset) 
     144     
     145    DO    l    = 1, llm 
     146      DO j=jj_begin-offset,jj_end+offset 
     147        DO i=ii_begin-offset,ii_end+offset 
     148          ij=(j-1)*iim+i 
     149          theta(ij,l) = theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) 
     150        ENDDO 
     151      ENDDO 
     152    ENDDO 
     153 
     154!$OMP BARRIER 
     155!$OMP MASTER 
     156      DEALLOCATE( p) 
     157!$OMP END MASTER 
     158!$OMP BARRIER     
     159     
     160  END SUBROUTINE compute_theta_rhodz2theta 
     161 
     162  SUBROUTINE compute_theta_rhodz2temperature(ps,theta_rhodz,temp,offset) 
     163  USE dimensions 
     164  USE geometry 
     165  USE metric 
     166  USE earth_const 
     167  USE pression_mod 
     168  USE exner_mod 
     169  IMPLICIT NONE 
     170    REAL(rstd),INTENT(IN) :: ps(iim*jjm) 
     171    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm) 
     172    REAL(rstd),INTENT(OUT) :: temp(iim*jjm,llm) 
     173    INTEGER,INTENT(IN) :: offset 
     174    INTEGER :: i,j,ij,l 
     175    REAL(rstd) :: p(iim*jjm,llm+1) 
     176    REAL(rstd) :: pk(iim*jjm,llm) 
     177    REAL(rstd) :: pks(iim*jjm) 
     178 
     179    CALL compute_pression(ps,p,offset) 
     180    CALL compute_exner(ps,p,pks,pk,offset) 
     181         
     182    DO    l    = 1, llm 
     183      DO j=jj_begin-offset,jj_end+offset 
     184        DO i=ii_begin-offset,ii_end+offset 
     185          ij=(j-1)*iim+i 
     186          temp(ij,l) = ( theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) ) * pk(ij,l) / cpp 
     187        ENDDO 
     188      ENDDO 
     189    ENDDO 
     190     
     191  END SUBROUTINE compute_theta_rhodz2temperature 
     192   
    58193END MODULE theta2theta_rhodz_mod 

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

    r12 r15  
    1616  USE transfert_mod 
    1717  USE metric 
    18   USE dissip_mod 
     18  USE dissip_gcm_mod 
    1919  USE ioipsl 
    2020  USE caldyn_gcm_mod 
    … …  
    2424  TYPE(t_field),POINTER :: f_phis(:) 
    2525  TYPE(t_field),POINTER :: f_theta(:) 
     26  TYPE(t_field),POINTER :: f_dtheta(:) 
    2627  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:) 
    2728  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:) 
    … …  
    5859  CALL getin('matsuno_period',matsuno_period) 
    5960  IF (scheme==leapfrog) matsuno_period=itaumax+1 
    60   
    61  
    6261 
    6362  CALL allocate_field(f_phis,field_t,type_real) 
    … …  
    7776  CALL allocate_field(f_dum2,field_u,type_real,llm) 
    7877 
     78  CALL allocate_field(f_theta,field_t,type_real,llm) 
     79  CALL allocate_field(f_dtheta,field_t,type_real,llm) 
     80 
    7981  CALL allocate_field(f_theta_rhodz,field_t,type_real,llm) 
    8082  CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm) 
    … …  
    8486  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm) 
    8587 
    86   CALL allocate_caldyn 
     88  CALL init_dissip(dt) 
     89  CALL init_caldyn(dt) 
    8790   
    8891!  CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u) 
    … …  
    9396    PRINT *,"It No :",It 
    9497 
    95     CALL caldyn(f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du) 
     98    CALL caldyn(it,f_phis,f_ps,f_theta_rhodz,f_u, f_dps, f_dtheta_rhodz, f_du) 
    9699 
    97100    IF (scheme==Euler) THEN 
    … …  
    102105      CALL  leapfrog_matsuno_scheme 
    103106    ELSE IF (scheme==adam_bashforth) THEN 
     107      CALL dissip(f_u,f_du,f_ps,f_theta_rhodz,f_dtheta_rhodz) 
    104108      CALL adam_bashforth_scheme 
    105109    ENDIF 

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

    r13 r15  
    1818  USE transfert_mod 
    1919  USE metric 
    20   USE dissip_mod 
     20  USE dissip_sw_mod 
    2121  USE ioipsl 
    2222 
    … …  
    8787      CALL  leapfrog_matsuno_scheme 
    8888    ELSE IF (scheme==adam_bashforth) THEN 
     89      CALL dissip(f_u,f_du)    
    8990      CALL adam_bashforth_scheme 
    9091    ENDIF 
    91  
    92     CALL dissip(f_u)    
    9392 
    9493  ENDDO 

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

    r12 r15  
    3131 
    3232    DO ind=1,ndomain 
     33      CALL swap_dimensions(ind) 
    3334      DO i=ii_begin,ii_end+1 
    3435        CALL request_add_point(ind,i,jj_begin-1,req_i1) 
    … …  
    6768    CALL create_request(field_u,req_e1) 
    6869    DO ind=1,ndomain 
     70      CALL swap_dimensions(ind) 
    6971      DO i=ii_begin,ii_end 
    7072        CALL request_add_point(ind,i,jj_begin-1,req_e1,rup) 

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

    r12 r15  
    402402                n=n+1 
    403403                CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n)) 
    404                 lon(n)=lon(n)*180/Pi+180 
     404                lon(n)=lon(n)*180/Pi 
    405405                lat(n)=lat(n)*180/Pi 
    406406                DO k=0,5 
    407407                  CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n)) 
    408408                  bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    409                   bounds_lon(k,n)=bounds_lon(k,n)*180/Pi+180 
     409                  bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
    410410                ENDDO 
    411411              ENDIF 
    … …  
    505505              n=n+1 
    506506              CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n)) 
    507               lon(n)=lon(n)*180/Pi+180 
     507              lon(n)=lon(n)*180/Pi 
     508 !             IF (lon(n)<0) lon(n)=lon(n)+360 
    508509              lat(n)=lat(n)*180/Pi 
    509510              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
    … …  
    513514              DO k=0,2 
    514515                bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    515                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi+180 
     516                bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     517!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
    516518              ENDDO 
    517519            ENDDO 
    … …  
    523525              n=n+1 
    524526              CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n)) 
    525               lon(n)=lon(n)*180/Pi+180 
     527              lon(n)=lon(n)*180/Pi 
     528!              IF (lon(n)<0) lon(n)=lon(n)+360 
    526529              lat(n)=lat(n)*180/Pi 
    527530              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n)) 
    … …  
    531534              DO k=0,2 
    532535                bounds_lat(k,n)=bounds_lat(k,n)*180/Pi 
    533                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi+180 
     536                bounds_lon(k,n)=bounds_lon(k,n)*180/Pi 
     537!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360 
    534538              ENDDO 
    535539            ENDDO 
Note: See TracChangeset for help on using the changeset viewer.