Changeset 286


Ignore:
Timestamp:
10/20/14 23:42:26 (10 years ago)
Author:
dubos
Message:

Partial etat0 cleanup (removed calls to xyz2lonlat)

Location:
codes/icosagcm/trunk/src
Files:
12 edited

Legend:

Unmodified
Added
Removed

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

    r266 r286  
    328328  INTEGER :: ij 
    329329  DO ij=ij_begin_ext,ij_end_ext 
    330      CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) 
    331      ulon(ij+u_right)=a*omega*cos(lat) 
     330     ulon(ij+u_right)=a*omega*cos(lat_e(ij+u_right)) 
    332331     ulat(ij+u_right)=0 
    333332 
    334      CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) 
    335      ulon(ij+u_lup)=a*omega*cos(lat) 
     333     ulon(ij+u_lup)=a*omega*cos(lat_e(ij+u_lup)) 
    336334     ulat(ij+u_lup)=0 
    337335 
    338      CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 
    339      ulon(ij+u_ldown)=a*omega*cos(lat) 
     336     ulon(ij+u_ldown)=a*omega*cos(lat_e(ij+u_ldown)) 
    340337     ulat(ij+u_ldown)=0 
    341338  END DO 

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

    r266 r286  
    256256          ij=(j-1)*iim+i 
    257257          IF (domain(ind)%own(i,j)) THEN 
    258 !            CALL xyz2lonlat(xyz_i(ij,:),lon,lat)  
    259 ! very bad inlining with intel compiler (>14 ?) 
    260              xyz(:)=xyz_i(ij,:) 
    261              xyz(:)=xyz(:)/sqrt(sum(xyz(:)**2)) 
    262              lat=asin(xyz(3)) 
    263              lon=atan2(xyz(2),xyz(1)) 
    264    
     258            lat=lat_i(ij)   
    265259            ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat)) 
    266260          END IF  

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

    r250 r286  
    540540 
    541541           nn = ij+shift_u 
    542            CALL xyz2lonlat(xyz_e(nn,:),lon,lat) 
    543542           zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm 
    544            CALL test2_schaer_mountain(lon,lat,p,z,zcoords,hybrid_eta, & 
     543           CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, & 
    545544                hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q) 
    546545!           u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:) 

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

    r266 r286  
    158158    REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) 
    159159 
    160     REAL(rstd) :: lon_i(iim*jjm) 
    161     REAL(rstd) :: lat_i(iim*jjm) 
    162160    REAL(rstd) :: temp_i(iim*jjm,llm) 
    163161    REAL(rstd) :: ulon_i(iim*jjm,llm) 
    164162    REAL(rstd) :: ulat_i(iim*jjm,llm) 
    165163 
    166     REAL(rstd) :: lon_e(3*iim*jjm) 
    167     REAL(rstd) :: lat_e(3*iim*jjm) 
    168164    REAL(rstd) :: ps_e(3*iim*jjm) 
    169165    REAL(rstd) :: mass_e(3*iim*jjm,llm) 
    … …  
    175171 
    176172    INTEGER :: l,i,j,ij 
    177  
    178     DO l=1,llm 
    179        DO j=jj_begin-1,jj_end+1 
    180           DO i=ii_begin-1,ii_end+1 
    181              ij=(j-1)*iim+i 
    182              CALL xyz2lonlat(xyz_i(ij,:)/radius,lon_i(ij),lat_i(ij)) 
    183              CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon_e(ij+u_right),lat_e(ij+u_right)) 
    184              CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon_e(ij+u_lup),lat_e(ij+u_lup)) 
    185              CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon_e(ij+u_ldown),lat_e(ij+u_ldown)) 
    186           END DO 
    187        END DO 
    188     END DO 
    189173 
    190174    SELECT CASE (TRIM(etat0_type)) 

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

    r186 r286  
    109109           DO i=ii_begin-1,ii_end+1 
    110110             ij=(j-1)*iim+i 
    111              CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
    112              ddsin=sin(lat)-sin(pi/20.) 
     111             ddsin=sin(lat_i(ij))-sin(pi/20.) 
    113112             theta(ij,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+thetarappel) 
    114113            ENDDO 
    … …  
    144143         ij=(j-1)*iim+i 
    145144 
    146          CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) 
     145         lat=lat_e(ij+u_right) 
    147146         IF (abs(sin(lat))<1.e-4) lat=1e-4 
    148147         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x 
    149148         fact(ij+u_right)=(1.-x)/(2.*omega*sin(lat)) 
    150149 
    151          CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) 
     150         lat=lat_e(ij+u_lup) 
    152151         IF (abs(sin(lat))<1.e-4) lat=1e-4 
    153152         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x 
    154153         fact(ij+u_lup)=(1.-x)/(2.*omega*sin(lat)) 
    155154 
    156          CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 
     155         lat=lat_e(ij+u_ldown) 
    157156         IF (abs(sin(lat))<1.e-4) lat=1e-4 
    158157         x=cos(lat) ; x=x*x ; x=x*x ; x=x*x 
    159158         fact(ij+u_ldown)=(1.-x)/(2.*omega*sin(lat)) 
    160  
    161 !         fact(ij+u_right)=-1 
    162 !         fact(ij+u_lup)=-1 
    163 !         fact(ij+u_ldown)=-1 
    164  
    165 !         CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 
    166 !         ps(ij)=lat 
    167159       ENDDO 
    168160     ENDDO 

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

    r271 r286  
    188188          DO i=ii_begin-1,ii_end+1 
    189189            n=(j-1)*iim+i 
    190             CALL xyz2lonlat(xyz_i(n,:),lon,lat) 
    191             CALL dist_lonlat(lon0,lat0,lon,lat,rr1)  ! GC distance from center  
     190            CALL dist_lonlat(lon0,lat0,lon_i(n),lat_i(n),rr1)  ! GC distance from center  
    192191            rr1 = radius*rr1 
    193192            IF ( rr1 .LT. R0 ) then  
    … …  
    212211          DO i=ii_begin-1,ii_end+1 
    213212            n=(j-1)*iim+i 
    214             CALL xyz2lonlat(xyz_i(n,:),lon,lat) 
    215             CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center  
     213            CALL dist_lonlat(lonc1,latc1,lon_i(n),lat_i(n),rr1)  ! GC distance from center  
    216214            rr1 = radius*rr1 
    217             CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center  
     215            CALL dist_lonlat(lonc2,latc2,lon_i(n),lat_i(n),rr2)  ! GC distance from center  
    218216            rr2 = radius*rr2 
    219217            dd1t1 = rr1/rt 
    … …  
    244242          DO i=ii_begin-1,ii_end+1 
    245243            n=(j-1)*iim+i 
    246             CALL xyz2lonlat(xyz_i(n,:),lon,lat) 
    247             CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center  
     244            CALL dist_lonlat(lonc1,latc1,lon_i(n),lat_i(n),rr1)  ! GC distance from center  
    248245            rr1 = radius*rr1 
    249             CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center  
     246            CALL dist_lonlat(lonc2,latc2,lon_i(n),lat_i(n),rr2)  ! GC distance from center  
    250247            rr2 = radius*rr2 
    251248            dd1t1 = rr1/rt 

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

    r186 r286  
    8686          DO i=ii_begin,ii_end 
    8787             ij=(j-1)*iim+i 
    88              CALL comp_all(xyz_i(ij,:), ps(ij),phis(ij),temp(ij,l), dummy1,dummy2) 
     88             CALL comp_all(lon_i(ij),lat_i(ij), ps(ij),phis(ij),temp(ij,l), dummy1,dummy2) 
    8989          END DO 
    9090       END DO 
    … …  
    100100          DO i=ii_begin,ii_end 
    101101             ij=(j-1)*iim+i 
    102              CALL comp_all(xyz_e(ij+u_right,:), dummy1,dummy2,dummy3, ulon(ij+u_right,l),ulat(ij+u_right,l)) 
    103              CALL comp_all(xyz_e(ij+u_lup,:), dummy1,dummy2,dummy3, ulon(ij+u_lup,l),ulat(ij+u_lup,l)) 
    104              CALL comp_all(xyz_e(ij+u_ldown,:), dummy1,dummy2,dummy3, ulon(ij+u_ldown,l),ulat(ij+u_ldown,l)) 
     102             CALL comp_all(lon_e(ij+u_right), lat_e(ij+u_right), & 
     103                  dummy1,dummy2,dummy3, ulon(ij+u_right,l),ulat(ij+u_right,l)) 
     104             CALL comp_all(lon_e(ij+u_lup), lat_e(ij+u_lup), & 
     105                  dummy1,dummy2,dummy3, ulon(ij+u_lup,l),ulat(ij+u_lup,l)) 
     106             CALL comp_all(lon_e(ij+u_ldown), lat_e(ij+u_ldown), & 
     107                  dummy1,dummy2,dummy3, ulon(ij+u_ldown,l),ulat(ij+u_ldown,l)) 
    105108          END DO 
    106109       END DO 
    … …  
    112115    CONTAINS 
    113116       
    114       SUBROUTINE comp_all(xyz, psj,phisj,tempj, ulonj,ulatj) 
     117      SUBROUTINE comp_all(lon,lat, psj,phisj,tempj, ulonj,ulatj) 
    115118        USE dcmip_initial_conditions_test_1_2_3 
    116         REAL(rstd), INTENT(IN) :: xyz(3) 
     119        REAL(rstd), INTENT(IN) :: lon, lat 
    117120        REAL(rstd), INTENT(OUT) :: psj,phisj,tempj,ulonj,ulatj 
    118         REAL :: lon,lat, dummy 
     121        REAL :: dummy 
    119122        dummy=0. 
    120         CALL xyz2lonlat(xyz,lon,lat) 
    121123        SELECT CASE (icase) 
    122124        CASE(0) 

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

    r186 r286  
    8383  REAL(rstd) :: Rd        ! gas constant of dry air, P=rho.Rd.T 
    8484  REAL(rstd) :: alpha, rm 
    85   REAL(rstd) :: lon,lat, C0, C1, GG 
     85  REAL(rstd) :: C0, C1, GG 
    8686  REAL(rstd) :: p0psk, pspsk,r,zz, thetab, thetap 
    8787  REAL(rstd) :: dummy, pp 
    … …  
    103103     DO i=ii_begin,ii_end 
    104104        ij=(j-1)*iim+i 
    105         CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
    106105 
    107106        IF(use_dcmip_routine) THEN 
    108            CALL test3_gravity_wave(lon,lat,pp,dummy,0, dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy) 
     107           CALL test3_gravity_wave(lon_i(ij),lat_i(ij),pp,dummy,0, dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy) 
    109108        ELSE 
    110            C1=C0*(cos(2*lat)-1) 
     109           C1=C0*(cos(2*lat_i(ij))-1) 
    111110            
    112111           !--- GROUND GEOPOTENTIAL 
    … …  
    120119            
    121120            
    122            r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)) 
     121           r=radius*acos(sin(latc)*sin(lat_i(ij))+cos(latc)*cos(lat_i(ij))*cos(lon_i(ij)-lonc)) 
    123122           s(ij)= d**2/(d**2+r**2) 
    124123        END IF 
    … …  
    134133           pp=0.5*(p(ij,l+1)+p(ij,l))  ! full-layer pressure 
    135134           IF(use_dcmip_routine) THEN 
    136               CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
    137               CALL test3_gravity_wave(lon,lat,pp,dummy,0,dummy,dummy,dummy,T(ij,l),dummy,dummy,dummy,dummy) 
     135              CALL test3_gravity_wave(lon_i(ij),lat_i(ij),pp,dummy,0, & 
     136                   dummy,dummy,dummy,T(ij,l),dummy,dummy,dummy,dummy) 
    138137           ELSE 
    139138              pspsk=(pp/ps(ij))**kappa 
    … …  
    162161           ij=(j-1)*iim+i 
    163162           IF(use_dcmip_routine) THEN 
    164               CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) 
    165               CALL test3_gravity_wave(lon,lat,pp,0.,0, & 
    166                    ulon(ij+u_right,l),ulat(ij+u_right,l),& 
     163              CALL test3_gravity_wave(lon_e(ij+u_right),lat_e(ij+u_right), & 
     164                   pp,0.,0, ulon(ij+u_right,l),ulat(ij+u_right,l),& 
    167165                   dummy,dummy,dummy,dummy,dummy,dummy) 
    168               CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) 
    169               CALL test3_gravity_wave(lon,lat,pp,0.,0,& 
    170                    ulon(ij+u_lup,l),ulat(ij+u_lup,l),& 
     166              CALL test3_gravity_wave(lon_e(ij+u_lup),lat_e(ij+u_lup), & 
     167                   pp,0.,0, ulon(ij+u_lup,l),ulat(ij+u_lup,l),& 
    171168                   dummy,dummy,dummy,dummy,dummy,dummy) 
    172               CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 
    173               CALL test3_gravity_wave(lon,lat,pp,0.,0,& 
    174                    ulon(ij+u_ldown,l),ulat(ij+u_ldown,l),& 
     169              CALL test3_gravity_wave(lon_e(ij+u_ldown),lat_e(ij+u_ldown), & 
     170                   pp,0.,0, ulon(ij+u_ldown,l),ulat(ij+u_ldown,l),& 
    175171                   dummy,dummy,dummy,dummy,dummy,dummy) 
    176172           ELSE 
    177               CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat) 
    178               ulon(ij+u_right,l)=u0*cos(lat) 
    179               ulat(ij+u_right,l)=0 
    180                
    181               CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat) 
    182               ulon(ij+u_lup,l)=u0*cos(lat) 
    183               ulat(ij+u_lup,l)=0 
    184                
    185               CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat) 
    186               ulon(ij+u_ldown,l)=u0*cos(lat) 
    187               ulat(ij+u_ldown,l)=0 
     173              ulon(ij+u_right,l) = u0*cos(lat_e(ij+u_right)) 
     174              ulat(ij+u_right,l) = 0 
     175              ulon(ij+u_lup,l)   = u0*cos(lat_e(ij+u_lup)) 
     176              ulat(ij+u_lup,l)   = 0 
     177              ulon(ij+u_ldown,l) = u0*cos(lat_e(ij+u_ldown)) 
     178              ulat(ij+u_ldown,l) = 0 
    188179           END IF 
    189180        ENDDO 

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

    r186 r286  
    121121          ij=(j-1)*iim+i 
    122122           
    123           CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat) 
     123          lon=lon_e(ij+u_right) ; lat=lat_e(ij+u_right) 
    124124          K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) 
    125125          r=radius*acos(K) 
    … …  
    127127          u(ij+u_right,l) = utot * sum(elon_e(ij+u_right,:) * ep_e(ij+u_right,:)) 
    128128 
    129  
    130           CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat) 
     129          lon=lon_e(ij+u_lup) ; lat=lat_e(ij+u_lup) 
    131130          K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) 
    132131          r=radius*acos(K) 
    … …  
    134133          u(ij+u_lup,l) = utot * sum(elon_e(ij+u_lup,:) * ep_e(ij+u_lup,:)) 
    135134 
    136           CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat) 
     135          lon=lon_e(ij+u_ldown) ; lat=lat_e(ij+u_ldown) 
    137136          K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc) 
    138137          r=radius*acos(K) 
    … …  
    151150         DO i=ii_begin-1,ii_end+1 
    152151           ij=(j-1)*iim+i 
    153            CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) 
    154              
     152           lat=lat_i(ij) 
    155153            Y(ij,l)=((-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5     & 
    156154                                + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega) 
    … …  
    168166       DO i=ii_begin,ii_end 
    169167         ij=(j-1)*iim+i 
    170          CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) 
     168         lat=lat_i(ij) 
    171169         phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sin(lat)**6 * (cos(lat)**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  & 
    172170                                           +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega ) 
    … …  
    191189           DO i=ii_begin,ii_end 
    192190             ij=(j-1)*iim+i 
    193              CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) 
     191             lon=lon_i(ij) ; lat=lat_i(ij) 
    194192             dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) &  
    195193                                        + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l)                    &  
    … …  
    220218           DO i=ii_begin,ii_end 
    221219             ij=(j-1)*iim+i 
    222              CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat) 
    223              q(ij,l,1)=q0*exp(-(lat/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2) 
     220             q(ij,l,1)=q0*exp(-(lat_i(ij)/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2) 
    224221           ENDDO 
    225222         ENDDO 

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

    r186 r286  
    66  TYPE(t_field),POINTER :: f_theta_eq(:) 
    77  TYPE(t_field),POINTER :: f_theta(:) 
    8   TYPE(t_field),POINTER :: f_clat(:) ! FIXME, duplication 
    98 
    109  REAL(rstd),ALLOCATABLE,SAVE :: knewt_t(:),kfrict(:) 
    … …  
    6564    REAL(rstd),POINTER :: u(:,:) 
    6665    REAL(rstd),POINTER :: q(:,:,:) 
    67     REAL(rstd),POINTER :: clat(:)  
    6866    REAL(rstd),POINTER :: theta_eq(:,:)  
    6967    REAL(rstd),POINTER :: theta(:,:)  
    … …  
    7876 
    7977       theta_eq=f_theta_eq(ind)  
    80        clat=f_clat(ind)  
    81        CALL compute_Teq(clat,theta_eq) ! FIXME : already done by Init_Teq 
     78       CALL compute_Teq(lat_i,theta_eq) ! FIXME : already done by Init_Teq 
    8279 
    8380       ps=f_ps(ind) 
    … …  
    111108       CALL allocate_field(f_theta,field_t,type_real,llm) 
    112109       CALL allocate_field(f_theta_eq,field_t,type_real,llm) 
    113        CALL allocate_field(f_clat,field_t,type_real) 
    114110       ALLOCATE(knewt_t(llm)); ALLOCATE( kfrict(llm))  
    115111 
    … …  
    147143          CALL swap_dimensions(ind) 
    148144          CALL swap_geometry(ind) 
    149           clat=f_clat(ind) 
    150145          theta_eq=f_theta_eq(ind) 
    151           CALL compute_Teq(clat,theta_eq) 
     146          CALL compute_Teq(lat_i,theta_eq) 
    152147       ENDDO 
    153148 
    … …  
    159154  END SUBROUTINE init_Teq 
    160155 
    161   SUBROUTINE compute_Teq(clat,theta_eq) 
     156  SUBROUTINE compute_Teq(lat,theta_eq) 
    162157    USE icosa 
    163158    USE disvert_mod 
    164159    IMPLICIT NONE   
    165     REAL(rstd),INTENT(OUT) :: clat(iim*jjm) 
     160    REAL(rstd),INTENT(IN) :: lat(iim*jjm) 
    166161    REAL(rstd),INTENT(OUT) :: theta_eq(iim*jjm,llm)  
    167162 
    168     REAL(rstd) :: lon, lat, r, zsig, ddsin, tetastrat, tetajl 
    169     REAL(rstd) :: slat(iim*jjm)  
     163    REAL(rstd) :: r, zsig, ddsin, tetastrat, tetajl 
    170164    INTEGER :: i,j,l,ij 
    171  
    172     DO j=jj_begin-1,jj_end+1 
    173        DO i=ii_begin-1,ii_end+1 
    174           ij=(j-1)*iim+i 
    175           CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
    176           clat(ij)=cos(lat)  
    177           slat(ij)=sin(lat)  
    178        ENDDO 
    179     ENDDO 
    180165 
    181166    DO l=1,llm 
    … …  
    185170          DO i=ii_begin-1,ii_end+1 
    186171             ij=(j-1)*iim+i 
    187              ddsin=slat(ij)  
     172             ddsin=SIN(lat(ij))  
    188173             tetajl=teta0-delt_y*ddsin*ddsin+eps*ddsin & 
    189174                  -delt_z*(1.-ddsin*ddsin)*log(zsig) 
    … …  
    244229       theta_eq=f_theta_eq(ind)  
    245230       theta=f_theta(ind)  
    246        clat=f_clat(ind)  
    247        CALL compute_heldsz(ps,theta_eq,clat, theta_rhodz,u, theta)  
     231       CALL compute_heldsz(ps,theta_eq,lat_i, theta_rhodz,u, theta)  
    248232    ENDDO 
    249233  END SUBROUTINE held_suarez 
    250234 
    251   SUBROUTINE compute_heldsz(ps,theta_eq,clat, theta_rhodz,u, theta)  
     235  SUBROUTINE compute_heldsz(ps,theta_eq,lat, theta_rhodz,u, theta)  
    252236    USE icosa 
    253237    USE theta2theta_rhodz_mod 
    … …  
    255239    REAL(rstd),INTENT(IN)    :: ps(iim*jjm)  
    256240    REAL(rstd),INTENT(IN)    :: theta_eq(iim*jjm,llm)  
    257     REAL(rstd),INTENT(IN)    :: clat(iim*jjm)  
     241    REAL(rstd),INTENT(IN)    :: lat(iim*jjm)  
    258242    REAL(rstd),INTENT(INOUT) :: theta_rhodz(iim*jjm,llm) 
    259243    REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm) 
    … …  
    268252             ij=(j-1)*iim+i 
    269253             theta(ij,l)=theta(ij,l) - dt*(theta(ij,l)-theta_eq(ij,l))* & 
    270                   (knewt_g+knewt_t(l)*clat(ij)**4 ) 
     254                  (knewt_g+knewt_t(l)*COS(lat(ij))**4 ) 
    271255          ENDDO 
    272256       ENDDO 

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

    r266 r286  
    33   
    44  TYPE t_geometry 
     5    TYPE(t_field),POINTER :: centroid(:) 
    56    TYPE(t_field),POINTER :: xyz_i(:) 
    6     TYPE(t_field),POINTER :: centroid(:) 
    77    TYPE(t_field),POINTER :: xyz_e(:) 
    88    TYPE(t_field),POINTER :: xyz_v(:) 
     9    TYPE(t_field),POINTER :: lon_i(:) 
     10    TYPE(t_field),POINTER :: lon_e(:) 
     11    TYPE(t_field),POINTER :: lat_i(:) 
     12    TYPE(t_field),POINTER :: lat_e(:) 
    913    TYPE(t_field),POINTER :: ep_e(:) 
    1014    TYPE(t_field),POINTER :: et_e(:) 
    … …  
    3034  REAL(rstd),POINTER :: Ai(:)          ! area of a cell 
    3135!$OMP THREADPRIVATE(Ai) 
     36  REAL(rstd),POINTER :: centroid(:,:)  ! coordinate of the centroid of the cell 
     37!$OMP THREADPRIVATE(centroid) 
    3238  REAL(rstd),POINTER :: xyz_i(:,:)     ! coordinate of the center of the cell (voronoi) 
    3339!$OMP THREADPRIVATE(xyz_i) 
    34   REAL(rstd),POINTER :: centroid(:,:)  ! coordinate of the centroid of the cell 
    35 !$OMP THREADPRIVATE(centroid) 
    3640  REAL(rstd),POINTER :: xyz_e(:,:)     ! coordinate of a wind point on the cell on a edge 
    3741!$OMP THREADPRIVATE(xyz_e) 
     42  REAL(rstd),POINTER :: xyz_v(:,:)     ! coordinate of a vertex (center of the dual mesh) 
     43!$OMP THREADPRIVATE(xyz_v) 
     44  REAL(rstd),POINTER :: lon_i(:)       ! longitude of the center of the cell (voronoi) 
     45!$OMP THREADPRIVATE(lon_i) 
     46  REAL(rstd),POINTER :: lon_e(:)       ! longitude of a wind point on the cell on a edge 
     47!$OMP THREADPRIVATE(lon_e) 
     48  REAL(rstd),POINTER :: lat_i(:)       ! latitude of the center of the cell (voronoi) 
     49!$OMP THREADPRIVATE(lat_i) 
     50  REAL(rstd),POINTER :: lat_e(:)       ! latitude of a wind point on the cell on a edge 
     51!$OMP THREADPRIVATE(lat_e) 
    3852  REAL(rstd),POINTER :: ep_e(:,:)      ! perpendicular unit vector of a edge (outsider) 
    3953!$OMP THREADPRIVATE(ep_e) 
    … …  
    4862  REAL(rstd),POINTER :: elat_e(:,:)    ! unit latitude vector on a wind point 
    4963!$OMP THREADPRIVATE(elat_e) 
    50   REAL(rstd),POINTER :: xyz_v(:,:)     ! coordinate of a vertex (center of the dual mesh) 
    51 !$OMP THREADPRIVATE(xyz_v) 
    5264  REAL(rstd),POINTER :: Av(:)          ! area of dual mesk cell 
    5365!$OMP THREADPRIVATE(Av) 
    … …  
    8496   
    8597    CALL allocate_field(geom%Ai,field_t,type_real,name='Ai') 
     98 
    8699    CALL allocate_field(geom%xyz_i,field_t,type_real,3) 
     100    CALL allocate_field(geom%lon_i,field_t,type_real) 
     101    CALL allocate_field(geom%lat_i,field_t,type_real) 
     102    CALL allocate_field(geom%elon_i,field_t,type_real,3) 
     103    CALL allocate_field(geom%elat_i,field_t,type_real,3) 
    87104    CALL allocate_field(geom%centroid,field_t,type_real,3) 
     105 
    88106    CALL allocate_field(geom%xyz_e,field_u,type_real,3) 
     107    CALL allocate_field(geom%lon_e,field_u,type_real) 
     108    CALL allocate_field(geom%lat_e,field_u,type_real) 
     109    CALL allocate_field(geom%elon_e,field_u,type_real,3) 
     110    CALL allocate_field(geom%elat_e,field_u,type_real,3) 
    89111    CALL allocate_field(geom%ep_e,field_u,type_real,3) 
    90112    CALL allocate_field(geom%et_e,field_u,type_real,3) 
    91     CALL allocate_field(geom%elon_i,field_t,type_real,3) 
    92     CALL allocate_field(geom%elat_i,field_t,type_real,3) 
    93     CALL allocate_field(geom%elon_e,field_u,type_real,3) 
    94     CALL allocate_field(geom%elat_e,field_u,type_real,3) 
     113 
    95114    CALL allocate_field(geom%xyz_v,field_z,type_real,3) 
    96115    CALL allocate_field(geom%de,field_u,type_real) 
    … …  
    119138    ep_e=geom%ep_e(ind) 
    120139    et_e=geom%et_e(ind) 
     140    lon_i=geom%lon_i(ind) 
     141    lat_i=geom%lat_i(ind) 
     142    lon_e=geom%lon_e(ind) 
     143    lat_e=geom%lat_e(ind) 
    121144    elon_i=geom%elon_i(ind) 
    122145    elat_i=geom%elat_i(ind) 
    … …  
    394417 
    395418          CALL xyz2lonlat(xyz_i(n,:),lon,lat) 
     419          lon_i(n)=lon 
     420          lat_i(n)=lat 
    396421          elon_i(n,1) = -sin(lon) 
    397422          elon_i(n,2) = cos(lon) 
    … …  
    403428           
    404429          CALL xyz2lonlat(xyz_e(n+u_right,:),lon,lat) 
     430          lon_e(n+u_right)=lon 
     431          lat_e(n+u_right)=lat 
    405432          elon_e(n+u_right,1) = -sin(lon) 
    406433          elon_e(n+u_right,2) = cos(lon) 
    … …  
    411438           
    412439          CALL xyz2lonlat(xyz_e(n+u_lup,:),lon,lat) 
     440          lon_e(n+u_lup)=lon 
     441          lat_e(n+u_lup)=lat 
    413442          elon_e(n+u_lup,1) = -sin(lon) 
    414443          elon_e(n+u_lup,2) = cos(lon) 
    … …  
    419448  
    420449          CALL xyz2lonlat(xyz_e(n+u_ldown,:),lon,lat) 
     450          lon_e(n+u_ldown)=lon 
     451          lat_e(n+u_ldown)=lat 
    421452          elon_e(n+u_ldown,1) = -sin(lon) 
    422453          elon_e(n+u_ldown,2) = cos(lon) 

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

    r214 r286  
    194194       CALL swap_geometry(ind) 
    195195       CALL pack_domain_2D(pack_info(ind), Ai, physics_inout%Ai) 
    196        CALL pack_lonlat(pack_info(ind)) 
     196       CALL pack_domain_2D(pack_info(ind), lon_i, physics_inout%lon) 
     197       CALL pack_domain_2D(pack_info(ind), lat_i, physics_inout%lat) 
    197198    END DO 
    198199  END SUBROUTINE init_pack_after 
    199  
    200   SUBROUTINE pack_lonlat(info) 
    201     USE icosa 
    202     IMPLICIT NONE 
    203     TYPE(t_pack_info) :: info 
    204     REAL(rstd) :: lon(iim*jjm), lat(iim*jjm) 
    205     INTEGER :: ij 
    206     DO ij=1,iim*jjm 
    207        CALL xyz2lonlat(xyz_i(ij,:),lon(ij),lat(ij))  
    208     ENDDO 
    209     CALL pack_domain_2D(info, lon, physics_inout%lon) 
    210     CALL pack_domain_2D(info, lat, physics_inout%lat) 
    211   END SUBROUTINE pack_lonlat 
    212200 
    213201!-------------------------------- Pack / Unpack 2D --------------------------- 
Note: See TracChangeset for help on using the changeset viewer.