subroutine testlambert() real :: gl, gb, x, y, lon0, lat0, y0, k, F, earth_radius, lat_stand1, GRIDWIDTH_M real :: deg2rad GRIDWIDTH_M = 2500.0 lon0 = 15.0 lat0 = 63.0 deg2rad = PI/180.0 earth_radius = 6371000.0 lat_stand1 = lat0 k = sin(PI/180.0*lat0) F = earth_radius*cos(PI/180.0*lat_stand1) * tan(PI/4 + PI/360.0*lat_stand1)**k/k y0 = F*tan(PI/4 - PI/360.0*lat0)**k gl = 15.0 gb = 63.0 call lb2lambert(x, y, gl, gb, lon0, y0, k, F) write(*,*) 'lon = ', gl, 'lat =', gb write(*,*) 'give lambert x = ', x, 'y =', y write(*,*) 'lambert i = ', (x)/GRIDWIDTH_M,'j =', y/GRIDWIDTH_M write(*,*) x = -892442.2 y = 1220678.0 call lambert2lb(x, y, gl, gb, lon0, y0, k, F) write(*,*) 'Lambert x = ', x, 'y =', y write(*,*) 'gives lon = ', gl, 'lat =', gb call lb2lambert(x, y, gl, gb, lon0, y0, k, F) write(*,*) 'and back to Lambert x = ', x, 'y =', y end subroutine testlambert